1

I am trying to round coordinates before querying a raster, to the precision of the raster, so that I can cache the results, and get the same values back even when the location is moved little. However, I find oddities in which pixel I access. In this case it's done through GDAL, but through Postgres it gives the same results.

gdalinfo ahn2_05m_dsm/r65dn1.tif:

Driver: GTiff/GeoTIFF
Files: ahn2_05m_dsm/r65dn1.tif
Size is 10000, 12500
Coordinate System is:
PROJCS["Amersfoort / RD New",
    GEOGCS["Amersfoort",
        DATUM["Amersfoort",
            SPHEROID["Bessel 1841",6377397.155,299.1528128,
                AUTHORITY["EPSG","7004"]],
            TOWGS84[565.4171,50.3319,465.5524,-0.398957388243134,0.343987817378283,-1.87740163998045,4.0725],
            AUTHORITY["EPSG","6289"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4289"]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.9999079],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    AUTHORITY["EPSG","28992"]]
Origin = (30000.000000000000000,393750.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   30000.000,  393750.000) (  3d35'10.68"E, 51d31'11.12"N)
Lower Left  (   30000.000,  387500.000) (  3d35'18.69"E, 51d27'48.94"N)
Upper Right (   35000.000,  393750.000) (  3d39'29.96"E, 51d31'15.04"N)
Lower Right (   35000.000,  387500.000) (  3d39'37.64"E, 51d27'52.85"N)
Center      (   32500.000,  390625.000) (  3d37'24.25"E, 51d29'32.01"N)
Band 1 Block=10000x1 Type=Float32, ColorInterp=Gray
  NoData Value=-3.4028234663852886e+38

You can see below that when the boundary for the pixel in the X direction is on .99|.00, however, in the Y direction it is .00|01 (actually it is .000000000|.0000000001).

$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30714.999999999 391055.00
Report:
  Location: (1429P,5390L)
  Band 1:
    Value: 5.90999984741211
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.00
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391054.99
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.00
Report:
  Location: (1430P,5390L)
  Band 1:
    Value: 5.71999979019165
$ gdallocationinfo -geoloc ahn2_05m_dsm/r65dn1.tif 30715.00 391055.01
Report:
  Location: (1430P,5389L)
  Band 1:
    Value: 6.1100001335144

Why is this, and is there a way to correct it?

2

This is happening because the geotransform sets the upper-left corner of the image. The 2D "pixel" interval is inclusive of the upper and left boundaries for each pixel and exclusive of the lower and right boundaries. The image interval isn't on the same coordinate system as the world coordinates (georeferencing), which makes for the "rounding difference" you are experiencing.

Let's say the origin is at x=100, y=200 and the pixel size is 1, -1 (see how the "Y" pixel size is negative? As you increase the line (y) you decrease the Y value.)

 (100, 200) (101, 200)       above this line is greater than 200
     v         v             and not part of the image.
     +---------+---------+   -----------------------------------
line |x v v v v|x v v v v|>
 |   |>        |>        |>
 v   |>        |>        |>
     |>        |>        |>  
     +---------+---------+ < (102, 199)
     |         |         |
     |         |         |
     |         |         |
     |         |         |
     +---------+---------+
           pixel ->

As for fixing this? The best option might be to mirror the image across a horizontal axis so that the new line_0, pixel_0 is at the lower-left corner of the geography you want to cover and update the geotransform accordingly. This would make the pixel/line coordinate system parallel to the X/Y world coordinate system.

1
  • Thanks! Had the feeling it had to have something to do with the flipped Y axis... – Tino Jul 20 '17 at 6:01

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.