3

My file is

Driver: GTiff/GeoTIFF
Files: 0001.TIF
Size is 8512, 8622
Coordinate System is:
PROJCS["WGS 84 / UTM zone 31N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32631"]]
Origin = (837903.261500865686685,5877808.477226150222123)
Pixel Size = (30.079329235631686,-30.079329235631686)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  837903.262, 5877808.477) (  8d 1'48.79"E, 52d56'35.71"N)
Lower Left  (  837903.262, 5618464.501) (  7d46'39.67"E, 50d37'11.83"N)
Upper Right ( 1093938.512, 5877808.477) ( 11d48'11.42"E, 52d43'20.22"N)
Lower Right ( 1093938.512, 5618464.501) ( 11d21'56.03"E, 50d24'59.35"N)
Center      (  965920.887, 5748136.489) (  9d44'38.52"E, 51d41'24.52"N)
Band 1 Block=8512x1 Type=UInt16, ColorInterp=Gray

After reprojecting it to EPSG:32632 by

gdalwarp -t_srs epsg:32632 0001.TIF 0001_.TIF

I got

Driver: GTiff/GeoTIFF
Files: 0001_.TIF
Size is 9192, 9293
Coordinate System is:
PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (413533.214507141499780,5866398.182513756677508)
Pixel Size = (29.999953551392952,-29.999953551392952)
Metadata:
  AREA_OR_POINT=Point
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  413533.215, 5866398.183) (  7d42'47.95"E, 52d56'24.88"N)
Lower Left  (  413533.215, 5587608.614) (  7d46'56.92"E, 50d26' 2.84"N)
Upper Right (  689292.788, 5866398.183) ( 11d48'55.18"E, 52d54'49.85"N)
Lower Right (  689292.788, 5587608.614) ( 11d39'51.06"E, 50d24'35.94"N)
Center      (  551413.001, 5727003.398) (  9d44'37.73"E, 51d41'29.87"N)
Band 1 Block=9192x1 Type=UInt16, ColorInterp=Gray

I calculated myself (using Java GeoTools) that

413533.215, 5608276.425 EPSG:32632 =

837903.262, 5618484.501 EPSG:32631

why I get such a big difference: 5587608.614 instead of 5608276.425? How does gdalwarp determine the pixel size and the extents?

EDIT:

The accepted answer shed light on my confusion. After I reprojected the lower right corner of the source raster instead of the lower left one, I got almost the same numbers:

x: 1093938.5119545625
y: 5618464.500556534
[668055.0297114132, 5587613.552738925]

x and y are the original coordinates in EPSG:32631

  • the valid area of one or the other of those projections probably fails to intersect the raster so you are likely to see all sorts of distortion and errors – Ian Turton Jun 30 '17 at 16:42
  • Probably, but I would like to know how gdalwarp does the calculations – Antonio Jun 30 '17 at 16:45
  • did you look at the code yet? – Ian Turton Jun 30 '17 at 16:46
  • I suppose there should be an algorithm or some considerations documented in a natural language... – Antonio Jun 30 '17 at 18:21
  • Gdaltransform gdal.org/gdaltransform.html gives 837903.261858582 5618464.50194627 as a result which proves that the algorithm has very little to do with your confusion. Please look at the image and accept the answer by @AndreJ. – user30184 Jun 30 '17 at 21:42
6

Just display your raster in a GIS software in UTM 32N:

enter image description here

You see that the raster is rotated in the new CRS, and the former Lower Left point is not the most southward point anymore.

Since gdalwarp always produces a rectangle parallel to the axes of the new projection, the new corner points are different from the UTM31N corner points.

| improve this answer | |

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.