Opened 16 years ago

Closed 9 years ago

#2117 closed defect (invalid)

Metadata Item GDALMD_AREA_OR_POINT is not used throughout.

Reported by: moellney Owned by: warmerdam
Priority: normal Milestone:
Component: GDAL_Raster Version: unspecified
Severity: normal Keywords:
Cc: Mateusz Łoskot

Description

It seems that the Metadata Item

GDALMD_AREA_OR_POINT

is not used correctly in GDAL formats.

Background:
For me it seems that GetGeoTransform(...) delivers the position of the "upper left corner" of the "raster data".

But the data may be points with extend (image) or without extend (DTM).

In the image case the geotransform information relates to the upper left corner of the upper left pixel. But if the raster data is interpreted as DTM values, then the points are just points: not the the center of an area (more the edge of an area)! GDAL does not interpret in GetGeoTransform(...), if the dataset should be an image or a DTM but in fact some raster formats give little hints. E.g. GeoTIFF and Esri ASCII Grid. These hints can be observed in MetaInformation GDALMD_AREA_OR_POINT.

For the DTM case it is important to interpret the geotransform right. Is it right that the Metainformation
GDALMD_AOP_POINT means: geotransform is position of upper left point in the dataset (and point has no area).
GDALMD_AOP_AREA means: geotransform is the position of the upper left corner of the area, which center is the upper right point in the dataset.

At least this would fit with the a definition I found for geotiff in http://www.remotesensing.org/geotiff/spec/geotiff2.5.html#2.5.2

But I found e.g. that in

USGSDEMDataset::LoadFromFile
...
adfGeoTransform[0] = (extent_min.x - dxdelta/2.0) / 3600.0;
...

the geotransform is calculated as if the points have an area to cover.

But in

GDALDataset *USGSDEMDataset::Open( GDALOpenInfo * poOpenInfo )
....
    poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
.....

we can read that the GDALMD_AREA_OR_POINT is set to point mode.

In

GDALDataset *AAIGDataset::Open( GDALOpenInfo * poOpenInfo )
...
    if ((i = CSLFindString( papszTokens, "xllcorner" )) >= 0 &&
        (j = CSLFindString( papszTokens, "yllcorner" )) >= 0 )
    {
        poDS->adfGeoTransform[0] = atof( papszTokens[i + 1] );
        poDS->adfGeoTransform[1] = dfCellDX;
        poDS->adfGeoTransform[2] = 0.0;
        poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
            + poDS->nRasterYSize * dfCellDY;
        poDS->adfGeoTransform[4] = 0.0;
        poDS->adfGeoTransform[5] = - dfCellDY;
    }
    else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 &&
             (j = CSLFindString( papszTokens, "yllcenter" )) >= 0 )
    {
        poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );

        poDS->adfGeoTransform[0] = atof(papszTokens[i + 1]) - 0.5 * dfCellDX;
        poDS->adfGeoTransform[1] = dfCellDX;
        poDS->adfGeoTransform[2] = 0.0;
        poDS->adfGeoTransform[3] = atof( papszTokens[j + 1] )
            - 0.5 * dfCellDY
            + poDS->nRasterYSize * dfCellDY;
        poDS->adfGeoTransform[4] = 0.0;
        poDS->adfGeoTransform[5] = - dfCellDY;
    }
...

For the case that we have XLLCORNER,YLLCORNER point mode is not set, but the geotransform is set equal to the metadata, no area is added. And for XLLCENTER,YLLCENTER an half cellsize area is added to geotranform...?

This is really confusing.

I add 4 small datasets in 4 different formats. The Dataset is always the same data (a 3x3 DTM values) that I took from geoengine.nima.mil for dted - so this is a DTM data that should be interpreted as points without area extend.

The upper left corner of the data is: x:10.983333 y: 51.0083333, the "distance" between the datapoints is: dx: 0.0166666 dy: -0.008333333

If you check the geotransform and the point/area metainfo flag you will see, that GTiff and IMG seems to work, but USGS DEM and ESRI ASCII deliver data that are not consistent:
ESRI: denotes Area but delivers the position of the point (no area extend added)
USGS DEM: denotes Point but delivers the position of the by area extended position of the point.

Isn't this all confusing.

I think a clear definition what POINT/AERA in the metainformation is missing for GDAL. Hey or I just did not find the place to look at.

Might all this be the reason for the fact that Google Earth has a half pixel shift in DTM??

See attached exmaples.

Attachments (4)

DTMs.zip (2.9 KB ) - added by moellney 16 years ago.
Four DTM formats of the same 3x3 points raster data.
WaterFlowingUpTheHill_China.kmz (722 bytes ) - added by moellney 16 years ago.
Faulty DTM in GoogleEarth (half a pixel shift for underlying SRTM 90m?)
WaterUpTheHill_Germany.kmz (717 bytes ) - added by moellney 16 years ago.
Also in Germany water can flow up the hills :-)
Referencing_Example.png (3.9 KB ) - added by moellney 16 years ago.
Illustration of anchor point for georeferencing

Download all attachments as: .zip

Change History (11)

by moellney, 16 years ago

Attachment: DTMs.zip added

Four DTM formats of the same 3x3 points raster data.

by moellney, 16 years ago

Faulty DTM in GoogleEarth (half a pixel shift for underlying SRTM 90m?)

by moellney, 16 years ago

Attachment: WaterUpTheHill_Germany.kmz added

Also in Germany water can flow up the hills :-)

comment:1 by warmerdam, 16 years ago

Component: defaultGDAL_Raster

I'm sorry, but I'm somewhat lost in this report. Is there a specific bug report you are making with regard to GDAL?

GDAL's geotransform is based on the assumption that pixels are areas. The GDALMD_AREA_OR_POINT metadata is only intended as a clue as to whether the values should be considered to have been collected over the whole area (area), or only at the areas center (point).

There has been much discussion about whether the GTiff driver ought to be adjusting the geotransform if the geotiff file is defined in pixelaspoint and basically my position is that pixel-is-point has no impact on georeferencing. It is just sensor interpretation. Lots of others disagree. If you wish this aspect of the GTiff driver to be your core bug report, then I believe you will find there is already a bug report open on it somewhere and that I'm unwilling to make any change at this time.

If you think that the semantics of GDALMD_AREA_OR_POINT need to be more clearly defined then lets narrow the focus of the bug report to that.

If you think there are errors in other file format drivers then please isolate one per bug report. For elevation based file formats like DTED I normally assume the origin values is the top left sampling point, and that the sample point in the GDAL model will be treated as applying over an area surrounding the sampling point, so I offset the origin by half a pixel to map to the GDAL model.

I suspect I need to write a GDAL FAQ on this topic so in the future I can just point to it rather getting incensed each time someone reports some aspect of this issue.

comment:2 by Mateusz Łoskot, 16 years ago

Cc: Mateusz Łoskot added

in reply to:  1 comment:3 by moellney, 16 years ago

Replying to warmerdam:

...munch...

GDAL's geotransform is based on the assumption that pixels are areas. The GDALMD_AREA_OR_POINT metadata is only intended as a clue as to whether the values should be considered to have been collected over the whole area (area), or only at the areas center (point).

Ok. I did not found this definition in the GDAL docs.

There has been much discussion about whether the GTiff driver ought to be adjusting the geotransform if the geotiff file is defined in pixelaspoint and basically my position is that pixel-is-point has no impact on georeferencing. It is just sensor interpretation. Lots of others disagree.

I think I agree with you. In GDAL we deal with rastered information. That means for me: we have data at positions (without extend) that are are 2-dim gridded.
To georeference this GDAL uses an affine transformation and needs an anchor point in relation of the data.

If I understand right this anchor point is the upper left corner of the upper left data point in the grid - however the upper left corner is defined by the extend of the area around each data point (please have a look at the attached image for this and correct me, if this is wrong).

...munch...

If you think that the semantics of GDALMD_AREA_OR_POINT need to be more clearly defined then lets narrow the focus of the bug report to that.

I did not find a place where is explicitly stated, that GDALMD_AREA_OR_POINT stands in no relation to the GetGeoTransform(...). I had the impression that it was related....

If you think there are errors in other file format drivers then please isolate one per bug report. For elevation based file formats like DTED I normally assume the origin values is the top left sampling point, and that the sample point in the GDAL model will be treated as applying over an area surrounding the sampling point, so I offset the origin by half a pixel to map to the GDAL model.

This is what I hoped to be. However the four DTM examples that I attached to this bug report (originating from geoengine.nima.mil) gave me two different georeferencing values back. These values could be reinterpreted with the GDALMD_AREA_OR_POINT information so that one georeferencing fits, but then a good on does not fit anymore.

I will fill bug reports for the formats that are do not create correct georeferencing values for this example ( I hope the examples are generated right).

I suspect I need to write a GDAL FAQ on this topic so in the future I can just point to it rather getting incensed each time someone reports some aspect of this issue.

Thanks for your clarification.

Michael

by moellney, 16 years ago

Attachment: Referencing_Example.png added

Illustration of anchor point for georeferencing

comment:4 by rayg, 16 years ago

I recommend dropping the metadata item and always using a pixels-as-areas protocol. There are only 12 refs to it in the entire lib, all of them at the driver level.

Rationale:

a) Drivers that encounter pixel-as-point data when reading can just quietly do an origin shift, and reverse shift if writing to pixel-as-point formats.

b) Would match Frank's existing intepretation, which is that the metadata didn't have any geospatial meaning anyway.

c) No new docs need be written; the existing dataset docs already assume pixels-as-areas.

in reply to:  4 comment:5 by moellney, 16 years ago

Replying to rayg:

I recommend dropping the metadata item and always using a pixels-as-areas protocol. There are only 12 refs to it in the entire lib, all of them at the driver level.

Rationale:

a) Drivers that encounter pixel-as-point data when reading can just quietly do an origin shift, and reverse shift if writing to pixel-as-point formats.

b) Would match Frank's existing intepretation, which is that the metadata didn't have any geospatial meaning anyway.

c) No new docs need be written; the existing dataset docs already assume pixels-as-areas.

I would like to support to (at least) rename the meta data GDALMD_AREA_OR_POINT to state more the difference between georeferencing and datacollection. Maybe GDALMD_AREA_OR_POINT_COLLECTED_DATA would be more distinguishable.

I got really confused (by probably wrong georeferenced files - half pixel shift problems) that could by fixed my misinterpreting the meta data GDALMD_AREA_OR_POINT. It looked to me as an exclusionary rule from the point_is_referenced_as_area default. And this came from the fact, that in the drivers the meta data GDALMD_AREA_OR_POINT is set in relation to the georeferencing coding.... as you can see in the original post for ARC ASCII grid.

in reply to:  1 comment:6 by moellney, 16 years ago

Replying to warmerdam: (...munch...)

GDAL's geotransform is based on the assumption that pixels are areas. The GDALMD_AREA_OR_POINT metadata is only intended as a clue as to whether the values should be considered to have been collected over the whole area (area), or only at the areas center (point).

Do you interpret the (x|y)llcorner or (x|y)llcenter as a hint if the data was collected area wide or only at the point.

I understand this tags as georeferencing tags as you do, too, when you shift the GDAL anchor according to the tags.

I could not find any definition of the ARC ASCII grid format, that supports the co-interpretation, that the data was collected in different ways regarding this tag:

...
    else if ((i = CSLFindString( papszTokens, "xllcenter" )) >= 0 &&
             (j = CSLFindString( papszTokens, "yllcenter" )) >= 0 )
    {
        poDS->SetMetadataItem( GDALMD_AREA_OR_POINT, GDALMD_AOP_POINT );
...

comment:7 by Even Rouault, 9 years ago

Resolution: invalid
Status: newclosed

Closing that ticket. This topic is always a hot one without a clear outcome...

Note: See TracTickets for help on using tickets.