search

opengeospatial / ogcapi-coverages

OGC API - Coverages draft specification
https://ogcapi.ogc.org/coverages
Apache License 2.0
22 stars 13 forks source link

Deal correctly with GeoTIFF Pixel-is-point vs Pixel-is-area #92

Open jratike80 opened 4 years ago

jratike80 commented 4 years ago

The difference between Pixel-is-point and Pixel-is-area in GeoTIFF tends to be confusing. I open this ticket for making my versions of the ASCII art that is used in the standard http://docs.opengeospatial.org/is/19-008r4/19-008r4.html and for giving some examples that can be re-produced with the provided test image and gdalinfo, gdalwarp, and listgeo commands.

GDAL and listgeo examples are given because GDAL is rather widely used and because it has some internal logic when it makes conversions between AREA and POINT variants. That logic is meant to provide backward interoperability after implementing the RFC https://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint and to hide some of the complexity from the average user. On the other hand, the same logic can make it difficult to understand what really happens. Fortunately listgeo and gdalinfo reveals all the secrets.

This is one of the ASCII illustrations from the standard

(0,0) (1,0)
      x---x---x-> I
      | * | * |
      x---x---x       Standard (PixelIsArea) TIFF Raster space R,
      | (1,1) (2,1)   showing the locations (x) and areas (*) of several pixels.
      v
      J
The location of the area pixel is referenced at the upper left of the pixel (pixelCorner convention).

These are my versions:

Case 1: Pixel is area. The ModelTiepointTags in the GeoTIFF are set to locations (0,0) and (5,5) in the raster space, and GTRasterTypeGeoKey is set to RasterPixelIsArea. Notice that the count of pixels is zero based and the lower right tiepoint in the 5 by 5 pixel image is actually the top left corner of the pixel at sixth row and column. It is OK, the tiepoint does not even need to be on the area that is covered by the tiff image and it does not need to match with pixel boundaries or centre.

pixel-is-area

Case 2: Pixel is point. The ModelTiepointTags in the GeoTIFF are set to the same locations (0,0) and (5,5) than in case 1 but GTRasterTypeGeoKey is set to RasterPixelIsPoint. It may feel more odd that the lower left tiepoint is now outside the surface that the raster covers but it is OK.

pixel-is-point

The side-by side presentation of these images as defined in the text of the standard

If a point-pixel image were to be displayed on a display device with pixel cells having the same size as the raster spacing, then the upper-left corner of the displayed image would be located in raster space at (-0.5, -0.5).

pixel-is-point-vs-pixel-is-area

Notice that the areas that cover all the pixels by the extreme boundaries are different.

The GDAL case and some examples

Traditionally RasterPixelIsArea vs RasterPixelIsPoint had no relevance to GDAL but that was changed at around the end of year 2010. But introducing half-a-pixel shifts to the extents of images by the value of the RasterTypeGeoKey is also confusing. Therefore the GDAL utilities like gdalinfo and gdal_translate are handling the RasterTypeGeoKey and ModelTiepointTags hand in hand.

The first example is a pixel-is-point GeoTIFF that has four tiepoints. Listgeo shows the tags

   Tagged_Information:
      ModelTiepointTag (8,3):
         0                 0                 0
         -180              90                0
         1024              0                 0
         -90               90                0
         0                 1024              0
         -180              0                 0
         1024              1024              0
         -90               0                 0
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeGeographic
      GTRasterTypeGeoKey (Short,1): RasterPixelIsPoint

Gdalinfo shows this

GCP[  0]: Id=1, Info=
          (-0.5,-0.5) -> (-180,90,0)
GCP[  1]: Id=2, Info=
          (1023.5,-0.5) -> (-90,90,0)
GCP[  2]: Id=3, Info=
          (-0.5,1023.5) -> (-180,0,0)
GCP[  3]: Id=4, Info=
          (1023.5,1023.5) -> (-90,0,0)
Metadata:
  AREA_OR_POINT=Point

Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1024.0)
Upper Right ( 1024.0,    0.0)
Lower Right ( 1024.0, 1024.0)
Center      (  512.0,  512.0)

Notice that because GDAL is internally using the tiepoint-at-corner approach it is automatically shifting the Ground Control Points that it is using by half a pixel when it sees the metadata that it is calling with name AREA_OR_POINT=Point.

Gdalinfo does not show georeferenced bounds for an image that has just tiepoints but we can gdalwarp the image and run gdalinfo for the warped north-up image

gdalwarp pixel_is_point.tif  warped_pixel_is_point.tif

gdalinfo warped_pixel_is_point.tif

Origin = (-179.956054687500000,89.956054687500000)
Pixel Size = (0.087890625000000,-0.087890625000000)
Metadata:
  AREA_OR_POINT=Point
Corner Coordinates:
Upper Left  (-179.9560547,  89.9560547) (179d57'21.80"W, 89d57'21.80"N)
Lower Left  (-179.9560547,  -0.0439453) (179d57'21.80"W,  0d 2'38.20"S)
Upper Right ( -89.9560547,  89.9560547) ( 89d57'21.80"W, 89d57'21.80"N)
Lower Right ( -89.9560547,  -0.0439453) ( 89d57'21.80"W,  0d 2'38.20"S)
Center      (-134.9560547,  44.9560547) (134d57'21.80"W, 44d57'21.80"N)

In the next test we can make a copy of the pixel_is_point image and edit just the GTRasterTypeGeoKey value

gdal_edit -mo AREA_OR_POINT=Area pixel_is_area.tif --config GTIFF_POINT_GEO_IGNORE YES

Now GDAL converts the tiepoints into its own GCPs without shifts and gdalinfo shows

GCP[  0]: Id=1, Info=
          (0,0) -> (-180,90,0)
GCP[  1]: Id=2, Info=
          (1024,0) -> (-90,90,0)
GCP[  2]: Id=3, Info=
          (0,1024) -> (-180,0,0)
GCP[  3]: Id=4, Info=
          (1024,1024) -> (-90,0,0)

Next use gdalwarp again and check gdalinfo

Origin = (-180.000000000000000,90.000000000000000)
Pixel Size = (0.087890625000000,-0.087890625000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
Lower Left  (-180.0000000,   0.0000000) (180d 0' 0.00"W,  0d 0' 0.01"N)
Upper Right ( -90.0000000,  90.0000000) ( 90d 0' 0.00"W, 90d 0' 0.00"N)
Lower Right ( -90.0000000,   0.0000000) ( 90d 0' 0.00"W,  0d 0' 0.01"N)
Center      (-135.0000000,  45.0000000) (135d 0' 0.00"W, 45d 0' 0.00"N)

Third test shows what GDAL is doing quietly by default. Make a new copy of the pixel-is-point image and edit it by using the GDAL defaults

gdal_edit -mo AREA_OR_POINT=Area pixel_is_area2.tif

I add here only the interesting part of the listgeo report.

  Tagged_Information:
      ModelTiepointTag (8,3):
         -0.5              -0.5              0
         -180              90                0
         1023.5            -0.5              0
         -90               90                0
         -0.5              1023.5            0
         -180              0                 0
         1023.5            1023.5            0
         -90               0                 0
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeGeographic
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea

Notice that by default GDAL is unwilling to slide the image. So what is does it to edit the GTRasterTypeGeoKey AND to slide the locations of the ModelTiepointTags by half a pixel. These two half-a-pixel slides compensate each other and therefore the pixel-is-area image that is produced this way overlaps exactly the original pixel-is-point image in QGIS, for example.

I guess that the automatic adjustment of georeferencing that GDAL is doing according to AREA_OR_POINT metadata can make people to believe that GDAL still does not care.

Test data and how to test:

pixel_vs_area.zip

jerstlouis commented 4 years ago

In this 256x256 pixel_is_point GeoTIFF tile:

https://maps.ecere.com/ogcapi/collections/SRTM_ViewFinderPanorama/tiles/GNOSISGlobalGrid/0/0/0

The upper-left corner is intended to represent the elevation at exactly lat: 90, lon -180; and the lower-right corner is intended to represent the elevation at exactly lat: -0.3515625, lon: -90.3515625

The TIE points right now are set to those coordinates in geographic space, and these in raster grid space:

1:   0.5,   0.5
2: 255.5,   0.5
3: 255.5,255.5
4:   0.5, 255.5

and gdalinfo reports:

Data axis to CRS axis mapping: 2,1
GCP[  0]: Id=1, Info=
          (0,0) -> (-180,90,0)
GCP[  1]: Id=2, Info=
          (255,0) -> (-90.3515625,90,0)
GCP[  2]: Id=3, Info=
          (0,255) -> (-180,0.3515625,0)
GCP[  3]: Id=4, Info=
          (255,255) -> (-90.3515625,0.3515625,0)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  256.0)
Upper Right (  256.0,    0.0)
Lower Right (  256.0,  256.0)
Center      (  128.0,  128.0)

The QGIS images show exactly side by side with https://maps.ecere.com/ogcapi/collections/SRTM_ViewFinderPanorama/tiles/GNOSISGlobalGrid/0/0/1

I am wondering whether I got this right, at least for pixel is point?

Thank you for all the clarifications!

joanma747 commented 4 years ago

My understanding of "pixel is point" is that the value (of elevation) is at the CENTER of each cell. This means that to have an elevation point at the north and south pole you need the cell to start a bit upper of 90 deg and end a bit lower of 90 deg.

jerstlouis commented 4 years ago

@joanma747 in my understanding, a measured value at an exact point does not cover an area (that's why this is Pixel Is POINT!), so I find this all extremely confusing. It is not necessarily the average of the values over an area... it is a sample at a precise point.

So this whole sliding by 0.5 really over-complicates things. Yes, when you want to render a map of a coverage you need to fill pixel occupying an area, so you could then assume that the sample point is the best approximation for an area equal to half a grid unit around the sample point, but this is only needed at the time that you "render" the data onto a map.

For example I may want to render this with dots of only 5x5 pixels regardless of how close I zoom, rather than square blocks occupying a space.

But now it seems that this 0.5 thing already has to be baked into the GeoTIFF, complicating even the simple data transfer and geo-referencing of those point samples, before even trying to render a map.

jratike80 commented 4 years ago

I started to think that GDAL is making somewhere a whole pixel error. Either my thinking is wrong or GDAL is sliding the image into wrong direction when the AREA_OR_POINT metadata is edited. I suggest to find some other problem to work with until I have talked with GDAL developers.

jerstlouis commented 4 years ago

@jratike80 My argument is that "sliding pixels" is completely wrong.

The user of the data simply needs to know whether a particular value represents an area or a point, and treat the value accordingly (e.g. when rendering a map). That is the only time that a half-grid-unit distance should be taken into account, if that is indeed how the tool wishes to display the data.

chris-little commented 4 years ago

@jerstlouis @joanma747 Pixel is Point and Pixel is Area are ideal cases. A Pixel is Point may be a 'point' measurement but the measuring process may be representative of the area. Remote sensing imagery is measuring an area, because of the optics and the finite size of the CCD sensors. It may not be an average either, but a max or min, because of signal thresholding, or some other 'cell-method' value.

Thinking about 'sliding' is not helpful either. Meteorologists think of two simultaneous grids, displaced by a half cell length. And we routinely use positions on the mid-points of the sides too. The various combinations were listed by Arakawa.

I have no problems with a half cell-length north of the North Pole or south of the South Pole., but have had loadsa problems with getting it right in software! Especially when we have moved the Equator to go through Paris France, and the North Pole is in the Pacific.

This doesn't help does it? ;-)

jerstlouis commented 4 years ago

@jratike80 @rouault Well a wrong direction shift would surely explain some of these things :)

jerstlouis commented 4 years ago

@chris-little I would argue that if it does represent an area (regardless of whether it is an average, min or max or other, due to finite size of CCD sensors or others), then Pixel is Area should be set.

In the case of Pixel is Point, the only assumption should be that this measurement is most exact at the exact geo-referenced point.

jratike80 commented 4 years ago

@jratike80 My argument is that "sliding pixels" is completely wrong.

The user of the data simply needs to know whether a particular value represents an area or a point, and treat the value accordingly (e.g. when rendering a map). That is the only time that a half-grid-unit distance should be taken into account, if that is indeed how the tool wishes to display the data.

You must play by the rules that you have given. GeoTIFF supports point vs area but the mechanism is connected to georeferencing. IMHO it does make some sense and is not completely wrong. Formats like png does not have such support and you must find your own way to inform your users about how to interpret the data. The strength of this approach is that it is easy to do it completely right.

Perhaps you like the GeoPackage gridded data extension http://docs.opengeospatial.org/is/17-066r1/17-066r1.html where pixels are stored raw and different metadata table stores the meaning with three alternatives "grid-value-is-center", "grid-value-is-area", or "grid-value-is-corner". It feels like "grid-value-is-corner" is what you would like to see in GeoTIFF. GeoTIFF is an OGC standard as well and perhaps the standard could be enhanced to match with GeoPackage here

http://www.opengis.net/spec/GeoTIFF/1.1/req/GTRasterTypeGeoKey.value The GTRasterTypeGeoKey value SHALL be:

Recommendation: the use of 0 (undefined) or 32767 (user-defined) is not recommended the only values

ghobona commented 3 years ago

From GML 3.2.2 "When a grid point is used to represent a sample space (e.g. image pixel), the grid point represents the center of the sample space (see ISO 19123:2005, 8.2.2)."

Note that the extract references ISO 19123 - Geographic information — Schema for coverage geometry and functions

jerstlouis commented 3 years ago

@jratike80

Perhaps you like the GeoPackage gridded data extension http://docs.opengeospatial.org/is/17-066r1/17-066r1.html where pixels are stored raw and different metadata table stores the meaning with three alternatives "grid-value-is-center", "grid-value-is-area", or "grid-value-is-corner". It feels like "grid-value-is-corner" is what you would like to see in GeoTIFF. GeoTIFF is an OGC standard as well and perhaps the standard could be enhanced to match with GeoPackage here

The GeoTIFF approach seems what is most correct to me (except for the fact that it also uses those 'offsets' which I find really confusing, considering that it already supports the Area vs Point distinction):

1 to indicate that the Raster type is PixelIsArea; or
2 to indicate that the Raster type is PixelIsPoint; or

assuming that PixelIsPoint means the GeoPackage tiled gridded coverage extension's grid-value-is-corner.

I believe grid-value-is-center is an offset in diguise which should generally be avoided, but in the case of the GeoPackage coverage tiles, it offers (maybe too much?) flexibility to select different samples when generating those tiles -- so it could be seen as PixelIsPoint in both cases, for which the coverage of the tiles has a different extent, and the distinction prescribes how to compute that extent.

If strictly applying the subset rules being proposed in http://github.com/opengeospatial/ogcapi-coverages/issues/95, typical TileMatrixSets based on 256x256 tiles would generate coverage tiles with 257x257 samples (including all 4 corner points) for PixelIsPoint, and 256x256 samples for PixelIsArea (with the exact same tile extents).

jerstlouis commented 2 years ago

SWG 2022-02-16: We should define a GeoTIFF encoding conformance class in OGC API - Coverages, and in there clarify that PixelIsArea / PixelIsPoint flags should be set accordingly, and that the data being returned from the subsetting operation should also accurately reflect the nature of the space occupied by the cells.

The subsetting rules should respect the "span" of the cells, and the a direct position (and the value at that position) should be included in the result set if the direct position is within the subset while considering the "span" of the cell for that value:

Note that if a GML response is negotiated, including application schemas that reference separate GeoTIFF or JPEG2000 rangesets (e.g. multi-part response), an "offset" vector is used from the CIS 1.0 logical model, but there is no way of knowing whether the data from a GML gmlcov:GridCoverage, gmlcov:RectifiedGridCoverage, gmlcov:ReferenceableGridCoverage cells represent point or area cell spaces. The valueIsPoint or valueIsArea can then only be indicated in e.g. the GeoTIFF.

In https://github.com/opengeospatial/coverage-implementation-schema/issues/6 we clarify how the distinction between point and area can be specified in CIS 1.1 GeneralGridCoverage.

prof-pguth commented 2 years ago

I have used GDAL and GRASS to create slope maps from a Pixel-is-point DEM. Both create a correct result, but shift the grid to Pixel-is-area and do not set the flag in Geotiff tag 1025. That tag is strongly encouraged in the Geotiff specs, which as I read them only allows two values, but does not specify one as the default.

https://github.com/OSGeo/grass/discussions/2261 has a long thread with the particulars in this case.

Can GDAL always set that flag in all Geotiffs it creates? That would remove any ambiguity. I agree with the results, just not the lack of specific details on how to interpret the shifted coordinates.

jratike80 commented 2 years ago

Please check your GDAL version. I believe that current GDAL version writes that tag always. The current version at the moment is 3.4.2.

I think that the question about "gdaldem slope" https://gdal.org/programs/gdaldem.html and the GeoTIFF tags should be asked from the gdal-dev mailing list. However, I made a quick test with GDAL 3.4.1, released 2021/12/27

gdaldem slope MNT_SRTM_EPSG32631.tif slope.tif
gdalinfo slope.tif
...
Metadata:
  AREA_OR_POINT=Area

listgeo slope.tif
...
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea

You have probably been reading the original GeoTIFF specification http://geotiff.maptools.org/spec/contents.html but GeoTIFF is now an OGC standard http://docs.opengeospatial.org/is/19-008r4/19-008r4.html with requirement 7.3

http://www.opengis.net/spec/GeoTIFF/1.1/req/GTRasterTypeGeoKey.value The GTRasterTypeGeoKey value SHALL be:

  • 0 to indicate that the Raster type is undefined or unknown; or
  • 1 to indicate that the Raster type is PixelIsArea; or
  • 2 to indicate that the Raster type is PixelIsPoint; or
  • 32767 to indicate that the Raster type is user-defined. Recommendation: the use of 0 (undefined) or 32767 (user-defined) is not recommended

If current GDAL version does not set GTRasterTypeGeoKey when it creates a new GeoTIFF then it feels like a bug but in that case you should make a re-producible test case for the develpers and write mail to gdal-dev list first.