Skip to content
This repository has been archived by the owner on Sep 1, 2022. It is now read-only.

Stereographic projection lat/lon bbox #504

Open
yuanho opened this issue Mar 28, 2016 · 5 comments
Open

Stereographic projection lat/lon bbox #504

yuanho opened this issue Mar 28, 2016 · 5 comments

Comments

@yuanho
Copy link
Member

yuanho commented Mar 28, 2016

In the GFS N Hemisphere dataset, the dataset bbox is ll: 21.65S 150.0W+ ur: 21.65S 30.00E which may not be wrong, but how can a client do any subset based on this value?

@cwardgar
Copy link
Contributor

The calculation for that lat-lon bounding box is done by ProjectionImpl.projToLatLonBB(ProjectionRect). It incorrectly assumes that you can calculate the lat-lon bounding box from the projection bounding box. I've got some code that demonstrates the problem:

  public static void main(String[] args) {
    // These projection parameters and values are from "Best GFS Northern Hemisphere 381km Time Series":
    // http:https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/N_Hemisphere_381km/Best.html
    Stereographic stereo = new Stereographic(90, 255, 0.9330127018922193);

    double[] x = new double[] {-12192.129, -11811.129, -11430.129, -11049.129, -10668.129, -10287.129, -9906.129,
            -9525.129, -9144.129, -8763.129, -8382.129, -8001.129, -7620.129, -7239.129, -6858.129, -6477.129,
            -6096.129, -5715.129, -5334.129, -4953.129, -4572.129, -4191.129, -3810.129, -3429.129, -3048.129,
            -2667.129, -2286.129, -1905.129, -1524.129, -1143.129, -762.129, -381.12903, -0.12901638, 380.87097,
            761.871, 1142.871, 1523.871, 1904.871, 2285.871, 2666.871, 3047.871, 3428.871, 3809.871, 4190.871,
            4571.871, 4952.871, 5333.871, 5714.871, 6095.871, 6476.871, 6857.871, 7238.871, 7619.871, 8000.871,
            8381.871, 8762.871, 9143.871, 9524.871, 9905.871, 10286.871, 10667.871, 11048.871, 11429.871, 11810.871,
            12191.871};
    double[] y = new double[] {-12192.129, -11811.129, -11430.129, -11049.129, -10668.129, -10287.129, -9906.129,
            -9525.129, -9144.129, -8763.129, -8382.129, -8001.129, -7620.129, -7239.129, -6858.129, -6477.129,
            -6096.129, -5715.129, -5334.129, -4953.129, -4572.129, -4191.129, -3810.129, -3429.129, -3048.129,
            -2667.129, -2286.129, -1905.129, -1524.129, -1143.129, -762.129, -381.12903, -0.12901638, 380.87097,
            761.871, 1142.871, 1523.871, 1904.871, 2285.871, 2666.871, 3047.871, 3428.871, 3809.871, 4190.871,
            4571.871, 4952.871, 5333.871, 5714.871, 6095.871, 6476.871, 6857.871, 7238.871, 7619.871, 8000.871,
            8381.871, 8762.871, 9143.871, 9524.871, 9905.871, 10286.871, 10667.871, 11048.871, 11429.871, 11810.871,
            12191.871};

    double minX = Double.MAX_VALUE;
    double minY = Double.MAX_VALUE;
    double maxX = -Double.MAX_VALUE;
    double maxY = -Double.MAX_VALUE;

    double minLon = Double.MAX_VALUE;
    double minLat = Double.MAX_VALUE;
    double maxLon = -Double.MAX_VALUE;
    double maxLat = -Double.MAX_VALUE;

    for (int i = 0; i < x.length; ++i) {
      minX = Math.min(minX, x[i]);
      minY = Math.min(minY, y[i]);
      maxX = Math.max(maxX, x[i]);
      maxY = Math.max(maxY, y[i]);

      ProjectionPoint projPoint = new ProjectionPointImpl(x[i], y[i]);
      LatLonPoint latLonPoint = stereo.projToLatLon(projPoint);
      System.out.printf("latLonPoint %s: %s%n", i + 1, latLonPoint);

      minLon = Math.min(minLon, latLonPoint.getLongitude());
      minLat = Math.min(minLat, latLonPoint.getLatitude());
      maxLon = Math.max(maxLon, latLonPoint.getLongitude());
      maxLat = Math.max(maxLat, latLonPoint.getLatitude());
    }

    ProjectionRect actualProjBB = new ProjectionRect(minX, minY, maxX, maxY);
    System.out.println("actualProjBB: " + actualProjBB);

    LatLonPoint minLatLon = new LatLonPointImpl(minLat, minLon);
    LatLonPoint maxLatLon = new LatLonPointImpl(maxLat, maxLon);
    LatLonRect actualLatLonBB = new LatLonRect(minLatLon, maxLatLon);
    System.out.println("actualLatLonBB: " + actualLatLonBB);

    LatLonRect wrongLatLonBB = stereo.projToLatLonBB(actualProjBB);
    System.out.println("wrongLatLonBB: " + wrongLatLonBB);
  }

The result:

latLonPoint 1: 20.8260S 150.0000W
latLonPoint 2: 19.1164S 150.0000W
latLonPoint 3: 17.3320S 150.0000W
latLonPoint 4: 15.4688S 150.0000W
latLonPoint 5: 13.5223S 150.0000W
latLonPoint 6: 11.4882S 150.0000W
latLonPoint 7: 9.3617S 150.0000W
latLonPoint 8: 7.1379S 150.0000W
latLonPoint 9: 4.8119S 150.0000W
latLonPoint 10: 2.3784S 150.0000W
latLonPoint 11: 0.1677N 150.0000W
latLonPoint 12: 2.8319N 150.0000W
latLonPoint 13: 5.6195N 150.0000W
latLonPoint 14: 8.5356N 150.0000W
latLonPoint 15: 11.5854N 150.0000W
latLonPoint 16: 14.7735N 150.0000W
latLonPoint 17: 18.1044N 150.0000W
latLonPoint 18: 21.5819N 150.0000W
latLonPoint 19: 25.2090N 150.0000W
latLonPoint 20: 28.9878N 150.0000W
latLonPoint 21: 32.9194N 150.0000W
latLonPoint 22: 37.0033N 150.0000W
latLonPoint 23: 41.2376N 150.0000W
latLonPoint 24: 45.6185N 150.0000W
latLonPoint 25: 50.1404N 150.0000W
latLonPoint 26: 54.7955N 150.0000W
latLonPoint 27: 59.5736N 150.0000W
latLonPoint 28: 64.4626N 150.0000W
latLonPoint 29: 69.4479N 150.0000W
latLonPoint 30: 74.5130N 150.0000W
latLonPoint 31: 79.6398N 150.0000W
latLonPoint 32: 84.8084N 150.0000W
latLonPoint 33: 89.9982N 150.0000W
latLonPoint 34: 84.8119N 30.0000E
latLonPoint 35: 79.6432N 30.0000E
latLonPoint 36: 74.5165N 30.0000E
latLonPoint 37: 69.4513N 30.0000E
latLonPoint 38: 64.4659N 30.0000E
latLonPoint 39: 59.5769N 30.0000E
latLonPoint 40: 54.7987N 30.0000E
latLonPoint 41: 50.1435N 30.0000E
latLonPoint 42: 45.6215N 30.0000E
latLonPoint 43: 41.2405N 30.0000E
latLonPoint 44: 37.0061N 30.0000E
latLonPoint 45: 32.9221N 30.0000E
latLonPoint 46: 28.9904N 30.0000E
latLonPoint 47: 25.2115N 30.0000E
latLonPoint 48: 21.5843N 30.0000E
latLonPoint 49: 18.1067N 30.0000E
latLonPoint 50: 14.7757N 30.0000E
latLonPoint 51: 11.5875N 30.0000E
latLonPoint 52: 8.5376N 30.0000E
latLonPoint 53: 5.6214N 30.0000E
latLonPoint 54: 2.8338N 30.0000E
latLonPoint 55: 0.1695N 30.0000E
latLonPoint 56: 2.3768S 30.0000E
latLonPoint 57: 4.8103S 30.0000E
latLonPoint 58: 7.1364S 30.0000E
latLonPoint 59: 9.3602S 30.0000E
latLonPoint 60: 11.4868S 30.0000E
latLonPoint 61: 13.5210S 30.0000E
latLonPoint 62: 15.4675S 30.0000E
latLonPoint 63: 17.3308S 30.0000E
latLonPoint 64: 19.1152S 30.0000E
latLonPoint 65: 20.8249S 30.0000E
actualProjBB: min: -12192.129 -12192.129 size: 24384.000 24384.000
actualLatLonBB:  ll: 20.8260S 150.0000W+ ur: 89.9982N 30.0000E
wrongLatLonBB:  ll: 20.8260S 150.0000W+ ur: 20.8249S 30.0000E

I'll turn this into a failing unit test. Thanks for the report, Yuan!

@ethanrd
Copy link
Member

ethanrd commented Mar 31, 2016

This dataset covers the entire Northern Hemisphere (plus some). All four corners of the grid are at 20.8S latitude with longitudes of 120E, 30E, 60W, and 150W.

The real problem is that a Latitude-Longitude Bounding Box just doesn't make sense when you are dealing with an area that includes one of the poles.

@dopplershift
Copy link
Member

Depends how you define "bounding box". To me that doesn't mean looking only at corners--that means looking at min/max extent. So if a pole is included, fine 90N is one of the ends. But by definition of bounds, unless that domain is actually a line in lon/lat space, you have a different value for the other end.

@cwardgar
Copy link
Contributor

My first code sample was wrong: you must include the lat/lon values of all grid cells in the bounding box calculation (the outer product of x and y). New code:

  public static void main(String[] args) {
    // These projection parameters and values are from "Best GFS Northern Hemisphere 381km Time Series":
    // http:https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/N_Hemisphere_381km/Best.html
    Stereographic stereo = new Stereographic(90, 255, 0.9330127018922193);

    double[] x = new double[] {-12192.129, -11811.129, -11430.129, -11049.129, -10668.129, -10287.129, -9906.129,
            -9525.129, -9144.129, -8763.129, -8382.129, -8001.129, -7620.129, -7239.129, -6858.129, -6477.129,
            -6096.129, -5715.129, -5334.129, -4953.129, -4572.129, -4191.129, -3810.129, -3429.129, -3048.129,
            -2667.129, -2286.129, -1905.129, -1524.129, -1143.129, -762.129, -381.12903, -0.12901638, 380.87097,
            761.871, 1142.871, 1523.871, 1904.871, 2285.871, 2666.871, 3047.871, 3428.871, 3809.871, 4190.871,
            4571.871, 4952.871, 5333.871, 5714.871, 6095.871, 6476.871, 6857.871, 7238.871, 7619.871, 8000.871,
            8381.871, 8762.871, 9143.871, 9524.871, 9905.871, 10286.871, 10667.871, 11048.871, 11429.871, 11810.871,
            12191.871};
    double[] y = new double[] {-12192.129, -11811.129, -11430.129, -11049.129, -10668.129, -10287.129, -9906.129,
            -9525.129, -9144.129, -8763.129, -8382.129, -8001.129, -7620.129, -7239.129, -6858.129, -6477.129,
            -6096.129, -5715.129, -5334.129, -4953.129, -4572.129, -4191.129, -3810.129, -3429.129, -3048.129,
            -2667.129, -2286.129, -1905.129, -1524.129, -1143.129, -762.129, -381.12903, -0.12901638, 380.87097,
            761.871, 1142.871, 1523.871, 1904.871, 2285.871, 2666.871, 3047.871, 3428.871, 3809.871, 4190.871,
            4571.871, 4952.871, 5333.871, 5714.871, 6095.871, 6476.871, 6857.871, 7238.871, 7619.871, 8000.871,
            8381.871, 8762.871, 9143.871, 9524.871, 9905.871, 10286.871, 10667.871, 11048.871, 11429.871, 11810.871,
            12191.871};

    double minLon = Double.MAX_VALUE;
    double minLat = Double.MAX_VALUE;
    double maxLon = -Double.MAX_VALUE;
    double maxLat = -Double.MAX_VALUE;

    for (int i = 0; i < x.length; ++i) {
      for (int j = 0; j < y.length; ++j) {
        ProjectionPoint projPoint = new ProjectionPointImpl(x[i], y[j]);
        LatLonPoint latLonPoint = stereo.projToLatLon(projPoint);

        minLon = Math.min(minLon, latLonPoint.getLongitude());
        minLat = Math.min(minLat, latLonPoint.getLatitude());
        maxLon = Math.max(maxLon, latLonPoint.getLongitude());
        maxLat = Math.max(maxLat, latLonPoint.getLatitude());
      }
    }

    ProjectionRect actualProjBB = new ProjectionRect(Doubles.min(x), Doubles.min(y), Doubles.max(x), Doubles.max(y));
    System.out.println("actualProjBB: " + actualProjBB);

    LatLonPoint minLatLon = new LatLonPointImpl(minLat, minLon);
    LatLonPoint maxLatLon = new LatLonPointImpl(maxLat, maxLon);
    LatLonRect actualLatLonBB = new LatLonRect(minLatLon, maxLatLon);
    System.out.println("actualLatLonBB: " + actualLatLonBB);

    LatLonRect wrongLatLonBB = stereo.projToLatLonBB(actualProjBB);
    System.out.println("wrongLatLonBB: " + wrongLatLonBB);
  }
actualProjBB: min: -12192.129 -12192.129 size: 24384.000 24384.000
actualLatLonBB:  ll: 20.8260S 179.9310W+ ur: 89.9982N 179.9323E
wrongLatLonBB:  ll: 20.8260S 150.0000W+ ur: 20.8249S 30.0000E

@JohnLCaron
Copy link
Collaborator

For the record, I think you are supposed to override ProjectionImpl.projToLatLonBB(ProjectionRect) to handle this. Note this is taking the projection coords, which should make sense, to latlon box, which as ethan says, may not.

The likely answer is you must subset in projection space, not latlon space.

Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

No branches or pull requests

5 participants