0

When I load a TanDEM GeoTIFF in ER Mapper, the coordinates of the cells are correct, i.e. the top left corner is the correct Lat/Long. When the same data are loaded in QGIS, the cells are all moved by 1/2 cell width (45m). The centre of the cell is the correct Lat/Long, not the corner.

I have checked the metadata for TanDEM, and each dataset should fit in an exact 1 degree square, so ER Mapper is correct, not QGIS.

I have to Warp the data in the same projection, with the correct XY max/min specified. When I reproject the data from WGS to UTM the error comes back, so I have to reproject again.

This looks like a bug in QGIS, as I can not find anywhere that allows me to specify cell centre or cell corner. This possibly needs to be an option to add unless someone can point me to an immediate solution. Am I missing an option somewhere?

2
  • The question is about this trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint, so pixel-is-point vs. pixel-is-area. Without having some sample data and documentation from the data provider it is hard to say which software is placing the image correctly. It is possible to change the interpretation of point vs. area but it is a bit tricky.
    – user30184
    Jan 20 '20 at 10:26
  • Thanks for the response and the link. I deal with satellite imagery a lot, and normally the coordinates represent the area (top corner) not the pixel centre. Here is part of the xml file that accompanies the TanDEM data, showing that the bounding box is an exact Lat/Long, while the pixel centre is 1/2 a pixel away: <westBoundingCoordinate>127.0000000 <southBoundingCoordinate>-19.0000000 <eastBoundingCoordinate>128.0000000 <northBoundingCoordinate>-18.0000000 <frameCoordsGeographic> <upperLeftLongitude>126.9995833 <upperLeftLatitude>-17.9995833 Jan 22 '20 at 6:48
2

TanDEM GeoTIFFs are of type "Pixel is point"

gdalinfo TDM1_DEM__30_N49E010_DEM.tif
Driver: GTiff/GeoTIFF
Files: TDM1_DEM__30_N49E010_DEM.tif
Size is 1201, 1201
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            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",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (9.999583333333334,50.000416666666666)
Pixel Size = (0.000833333333333,-0.000833333333333)
Metadata:
  AREA_OR_POINT=Point
  TIFFTAG_DATETIME=2014:12:25 18:30:34
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=W42/DLR
  TIFFTAG_XRESOLUTION=50
  TIFFTAG_YRESOLUTION=50
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   9.9995833,  50.0004167) (  9d59'58.50"E, 50d 0' 1.50"N)
Lower Left  (   9.9995833,  48.9995833) (  9d59'58.50"E, 48d59'58.50"N)
Upper Right (  11.0004167,  50.0004167) ( 11d 0' 1.50"E, 50d 0' 1.50"N)
Lower Right (  11.0004167,  48.9995833) ( 11d 0' 1.50"E, 48d59'58.50"N)
Center      (  10.5000000,  49.5000000) ( 10d30' 0.00"E, 49d30' 0.00"N)
Band 1 Block=1201x2 Type=Float32, ColorInterp=Gray
  Unit Type: metre

QGIS behaves in the same way than GDAL since the implementation of RFC 33 https://trac.osgeo.org/gdal/wiki/rfc33_gtiff_pixelispoint

Planned Changes

Interpretation of the raster space from the GeoTIFF tie points will be offset by half a pixel in the PixelIsPoint case in gdal/frmts/gtiff/geotiff.cpp. This will impact the formation of the geotransform and the formation of GCPs when there are multiple tie points. geotransmatrix conversion to geotransform will also be affected. Conversely if writing files with PixelIsPoint (as driven by the "AREA_OR_POINT" metadata item being set to "POINT") the written raster space coordinates would be offset by half a pixel. In trunk the above behavior may be disabled by setting the GTIFF_POINT_GEO_IGNORE configuration option to TRUE (it will default to FALSE). In GDAL 1.7 and 1.6 branch the same changes will be applied, except the GTIFF_POINT_GEO_IGNORE configuration option will default to TRUE.

If you want to fix the image to behave as you wish you can do it with gdal_edit https://gdal.org/programs/gdal_edit.html.

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

Now the georeferencing is different

Corner Coordinates:
Upper Left  (  10.0000000,  50.0000000) ( 10d 0' 0.00"E, 50d 0' 0.00"N)
Lower Left  (  10.0000000,  48.9991667) ( 10d 0' 0.00"E, 48d59'57.00"N)
Upper Right (  11.0008333,  50.0000000) ( 11d 0' 3.00"E, 50d 0' 0.00"N)
Lower Right (  11.0008333,  48.9991667) ( 11d 0' 3.00"E, 48d59'57.00"N)
Center      (  10.5004167,  49.4995833) ( 10d30' 1.50"E, 49d29'58.50"N)

Just keep in mind that fundamentally the image has now wrong georeferencing.

Edit: The image metadata XML contains both extents. The boundingBox is measused from the centers of the pixels while the frameCoordsGeographic make a box that is half a pixel wider in all directions. You are right that half a pixel shifts are tricky and with a DEM having big pixel size it may matter if the height of the terrain is measured from the centre vs. from the corner of a pixel. On the other hand, for analysis it may be more important to make the pixels to aline and thus prevent resampling. Just know your data and make reasonable decisions.

<spatialCoverage>
            <boundingBox>
            <westBoundingCoordinate>10.0000000</westBoundingCoordinate>
            <southBoundingCoordinate>49.0000000</southBoundingCoordinate>
            <eastBoundingCoordinate>11.0000000</eastBoundingCoordinate>
            <northBoundingCoordinate>50.0000000</northBoundingCoordinate>
            </boundingBox>
            <frameCoordsGeographic>
                <upperLeftLongitude>9.9995833</upperLeftLongitude>
                <upperLeftLatitude>50.0004167</upperLeftLatitude>
                <upperRightLongitude>11.0004167</upperRightLongitude>
                <upperRightLatitude>50.0004167</upperRightLatitude>
                <lowerRightLongitude>11.0004167</lowerRightLongitude>
                <lowerRightLatitude>48.9995833</lowerRightLatitude>
                <lowerLeftLongitude>9.9995833</lowerLeftLongitude>
                <lowerLeftLatitude>48.9995833</lowerLeftLatitude>
            </frameCoordsGeographic>
1
  • Thanks for the very helpful reply. The TanDEM data needs to fit in the bounding box given in the xml file, and that is exact Lat/Long. I can see that there are ramifications if changes are made to the algorithm. My workaround is OK, so my main point is that data users need to be really careful about where they think there data lies in space. I have verified the TanDEM data against other sources (ALOS, SRTM. GPS tracklogs and Landsat), so am confident everything lines up. Jan 22 '20 at 6:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.