search

opengeospatial / CRS-Gridded-Geodetic-data-eXchange-Format

Gridded Geodetic data eXchange Format
11 stars 3 forks source link

Dealing with NADCON5 geographic3dOffsets grids whose resolution for ellipsoidal height difference is != lat and lon #52

Open rouault opened 1 year ago

rouault commented 1 year ago

Just looking at NADCON5 grids, I've found that some of those grids have a different resolution for the ellipsoidal height difference than for the lat & lon differences.

For example https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/nad83_1992.nad83_2007.alaska/nadcon5.nad83_1992.nad83_2007.alaska.eht.trn.20160901.b has a 0.25 x 0.25 degree resolution (241 x 93 pixels )whereas https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/nad83_1992.nad83_2007.alaska/nadcon5.nad83_1992.nad83_2007.alaska.lat.trn.20160901.b has a 0.125 x 0.125 degree resolution (481 x 185 pixels) Luckily they have the same spatial extent though. Given the above, if using the raw data, one cannot create a single geographic3dOffsets grid.

One possibility would be for the grid producer to interpolate to the ellipsoidal height difference grid to upgrade it to 0.125 x 0.125 degree resolution. But, NADCON5 uses biquadratic interpolation, and, contrary to bilinear interpolation, I wouldn't expect the result of biquadratic interpolation on a biquadatic interpolated grid to be exactly the same as biquadric interpolation on the original grid.

Or should that be stored as 2 separate grids in a file, one with geographic2dOffsets and another one with a content_type=EllipsoidalHeightOffset (not currently available)

ccrook commented 1 year ago

Interesting question. The latter option would work (at least with the proof of concept implementation), ie two grids each in its own ggxfGroup, one with parameters latitudeOffset, longitudeOffset, the other with just ellipsoidalHeightOffset. These are evaluated separately, each with its own grid and each setting the not included parameters to zero, and then added.

This does raise a question for me though, which is that we haven't really documented how multiple GGXF groups are to be used. For deformation models the above approach is used, and as they include time functions these are calculated at the coordinate epoch and multiply displacements for the ggxfGroups before they are summed. In the proof of concept code the multiplier is 1.0 if a group doesn't include a time function. However that was all put in the POC code early on to demonstrate that it could correctly calculate the values, and therefore that the format was serving its purpose.

So far this hasn't been documented in the GGXF spec. It is defined in the deformation specification but that doesn't explicitly define how it is mapped to GGXF.

A very interesting test case Even, thanks.

ccrook commented 1 year ago

@rouault I was hoping to look at these with gdal but looks like it was added very recently as NOAA_B. I tried osgeo/gdal docker image, but that appears to be at 3.6. If you have time (!) please could you convert to another format suggest easy way to read! Thanks

rouault commented 1 year ago

it was added very recently as NOAA_B

yes, just a few days ago in master.

please could you convert to another format suggest easy way to read!

yes, as my favorite flavor of GGXF, aka GTG: us_noaa_nadcon5_nad83_1992_nad83_2007_alaska.tif.zip. It contains 2 grids: one for lat-lon offsets and one for ellipsoidal height offset

$ gdalinfo GTIFF_DIR:1:us_noaa_nadcon5_nad83_1992_nad83_2007_alaska.tif
Driver: GTiff/GeoTIFF
Files: us_noaa_nadcon5_nad83_1992_nad83_2007_alaska.tif
Size is 481, 185
Coordinate System is:
GEOGCRS["NAD83(HARN)",
    DATUM["NAD83 (High Accuracy Reference Network)",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4152]]
Data axis to CRS axis mapping: 2,1
Origin = (171.937500000000000,73.062500000000000)
Pixel Size = (0.125000000000000,-0.125000000000000)
Metadata:
  area_of_use=USA - Alaska
  AREA_OR_POINT=Point
  auxiliary_data=ellipsoidal_height_offset
  interpolation_method=biquadratic
  target_crs_epsg_code=4893
  TIFFTAG_COPYRIGHT=Derived from work by NOAA. Public Domain
  TIFFTAG_DATETIME=2016:09:01 00:00:00
  TIFFTAG_IMAGEDESCRIPTION=EPSG:4152 (NAD83(HARN)) to EPSG:4893 (NAD83(NSRS2007)). Converted from https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/nad83_1992.nad83_2007.alaska
  TYPE=HORIZONTAL_OFFSET
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=PIXEL
  PREDICTOR=3
Corner Coordinates:
Upper Left  ( 171.9375000,  73.0625000) (171d56'15.00"E, 73d 3'45.00"N)
Lower Left  ( 171.9375000,  49.9375000) (171d56'15.00"E, 49d56'15.00"N)
Upper Right (     232.062,      73.062) (232d 3'45.00"E, 73d 3'45.00"N)
Lower Right (     232.062,      49.938) (232d 3'45.00"E, 49d56'15.00"N)
Center      (     202.000,      61.500) (202d 0' 0.00"E, 61d30' 0.00"N)
Band 1 Block=256x256 Type=Float32, ColorInterp=Gray
  Description = latitude_offset
  Unit Type: arc-second
Band 2 Block=256x256 Type=Float32, ColorInterp=Undefined
  Description = longitude_offset
  Unit Type: arc-second
  Metadata:
    positive_value=east
$ gdalinfo GTIFF_DIR:2:us_noaa_nadcon5_nad83_1992_nad83_2007_alaska.tif
Driver: GTiff/GeoTIFF
Files: us_noaa_nadcon5_nad83_1992_nad83_2007_alaska.tif
Size is 241, 93
Coordinate System is:
GEOGCRS["NAD83(HARN)",
    DATUM["NAD83 (High Accuracy Reference Network)",
        ELLIPSOID["GRS 1980",6378137,298.257222101004,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4152]]
Data axis to CRS axis mapping: 2,1
Origin = (171.875000000000000,73.125000000000000)
Pixel Size = (0.250000000000000,-0.250000000000000)
Metadata:
  area_of_use=USA - Alaska
  AREA_OR_POINT=Point
  auxiliary_data=horizontal_offset
  interpolation_method=biquadratic
  target_crs_epsg_code=4893
  TIFFTAG_COPYRIGHT=Derived from work by NOAA. Public Domain
  TIFFTAG_DATETIME=2016:09:01 00:00:00
  TIFFTAG_IMAGEDESCRIPTION=EPSG:4152 (NAD83(HARN)) to EPSG:4893 (NAD83(NSRS2007)). Converted from https://geodesy.noaa.gov/pub/nadcon5/20160901release/Builds/nad83_1992.nad83_2007.alaska
  TYPE=ELLIPSOIDAL_HEIGHT_OFFSET
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
  PREDICTOR=3
Corner Coordinates:
Upper Left  ( 171.8750000,  73.1250000) (171d52'30.00"E, 73d 7'30.00"N)
Lower Left  ( 171.8750000,  49.8750000) (171d52'30.00"E, 49d52'30.00"N)
Upper Right (     232.125,      73.125) (232d 7'30.00"E, 73d 7'30.00"N)
Lower Right (     232.125,      49.875) (232d 7'30.00"E, 49d52'30.00"N)
Center      (     202.000,      61.500) (202d 0' 0.00"E, 61d30' 0.00"N)
Band 1 Block=241x93 Type=Float32, ColorInterp=Gray
  Description = ellipsoidal_height_offset
  Unit Type: metre

There are other NADCON5 grids where the grid with elliposidal height offset is actually more resoluted than the lat-lon offset one.

ccrook commented 1 year ago

Thanks Even - perfect. Took me a while to figure out just how recent it was after I saw the driver listed on the GDAL raster drivers page! Also thanks for your help with clarifying the relationship of the Dataset ReadAsArray and GetGeoTransform axis orders.

RogerLott commented 1 year ago

@ccrook In the geographic3dOffsets method supported by GGXF, there is no addition or combination of the three offsets evaluated for a point. Issue #53 is not needed to handle the question here. In the geographic3dOffsets method the offsets are applied through the three independent formulas defined by the method.

In GGXF, grids with different resolution have to be handled separately. One possibility is as you suggest, using different ggxfGroups.

However the practical way of treating them will probably be in three grids within one ggxfGroup. The different grid spacing comes out in the wash through the array of affine coefficients and node count for each grid:

!! Care - affine coeffs here are for illustration. The values for coefficients A1 and B1 are ficticious, and depending upon corner used as grid origin the signs for coefficients A3 and B2 might be transposed. The illustration is that the different node spacing is embedded in coefficients A3 and B2. Also note that GGXF requires longitudes and longitude offsets positive east, and current NADCON longitudes are positive west, so the direction reversal will need to be treated in the GGXF binary file creation.

ccrook commented 1 year ago

I think the structure would look more like the following. The horizontal and vertical grids have different node parameters (ie latitudeOffset, longitudeOffset for the horizontal and ellipsoialHeightOffset for the vertical), so they must be in different ggxfGroups. (If the horizontal grid is split into lat and lon grids then they would also be in separate ggxfGroups, however as they have the same grid spacing there is no need to do so)

ggxfGroups:
- ggxfGroupName: "horizontal_offset"
  interpolationMethod: biquadratic
  gridParameters:
  - latitudeOffset
  - longitudeOffset
  grids:
  - gridName: alaska-horizonal-offset
    affineCoeffs: [73.0, 0.0, -0.125, 172.0, 0.125, 0.0]
    iNodeCount: 481
    jNodeCount: 185
    dataSource:
       ....
- ggxfGroupName: "vertical_offset"
  interpolationMethod: biquadratic
  gridParameters:
  - ellipsoidalHeightOffset
  grids:
  - gridName: alaska-vertical-offset
    affineCoeffs: [73.0, 0.0, -0.25, 172.0, 0.25, 0.0]
    iNodeCount: 241
    jNodeCount: 93
    dataSource:
     ....