Monday 31 July 2017

geotiff tiff - Why doesn't gdalinfo report pixel size?


When I run gdalinfo (v2.2.2) on a geotiff that I created from a pdf using gdal_translate, I don't get the pixel size as shown in the docs. Does the pixel size have to be explicitly stored in the tif in order for gdalinfo to report it? If so, can gdal_translate be coerced to store the value when the tif is written?



>gdalinfo NM_Canada_Ojitos_20170216.tif
Driver: GTiff/GeoTIFF
Files: NM_Canada_Ojitos_20170216.tif
Size is 3412, 4350
Coordinate System is:
PROJCS["NAD83 / UTM zone 13N",
GEOGCS["NAD83",
DATUM["North_American_Datum_1983",
SPHEROID["GRS 1980",6378137,298.257222101,
AUTHORITY["EPSG","7019"]],

TOWGS84[0,0,0,0,0,0,0],
AUTHORITY["EPSG","6269"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4269"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-105],

PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","26913"]]
GeoTransform =
319569.9072199312, 4.06320686118055, -0.08151911420909382

4042818.402645736, -0.08151911420909382, -4.06320686118055
Metadata:
AREA_OR_POINT=Area
AUTHOR=USGS National Geospatial Technical Operations Center
CREATION_DATE=D:20170216103259Z
CREATOR=ESRI ArcSOC 10.0.2.3200
KEYWORDS=Topographic, Transportation, Hydrography, Orthoimage, U.S. National Grid, imageryBaseMapsEarthCover, Imagery and Base Ma
s, Geographic Names Information System
NEATLINE=POLYGON ((332137.201278231 4041102.99001007,331856.525171518 4027113.079869,320528.939318619 4027340.34242206,320809.615
25332 4041330.25256312,332137.201278231 4041102.99001007))

SUBJECT=This image map depicts geographic features on the surface of the earth. It was created to provide a representation of ac
essible geospatial data which is readily available to enhance the capability of Federal, State, and local emergency responders for
omeland security efforts. This image map is generated from selected National Map data holdings and other cartographic data.
TITLE=USGS 7.5-minute image map for Canada Ojitos, New Mexico
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 319569.907, 4042818.403) (107d 0'53.88"W, 36d30'49.40"N)
Lower Left ( 319215.299, 4025143.453) (107d 0'53.29"W, 36d21'15.87"N)

Upper Right ( 333433.569, 4042540.259) (106d51'36.58"W, 36d30'49.43"N)
Lower Right ( 333078.961, 4024865.310) (106d51'37.13"W, 36d21'15.87"N)
Center ( 326324.434, 4033841.856) (106d56'15.22"W, 36d26' 2.73"N)
Band 1 Block=3412x1 Type=Byte, ColorInterp=Red
Band 2 Block=3412x1 Type=Byte, ColorInterp=Green
Band 3 Block=3412x1 Type=Byte, ColorInterp=Blue

Answer



The issue is that your raster is rotated. gdalinfo intentionally only reports origin and pixel size for north up rasters, it will show the full geotransform for rotated rasters.


You can extract the pixel sizes from the geotransform with a little maths (see below) so I'm not sure why GDAL doesn't report the pixel size anyway. Maybe because it might be confusing as it will be different to the values shown in the geotransform.


You can make your raster north up by running gdalwarp on it and then gdalinfo will give you pixel size:



$ gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
Size is 3412, 4350
Coordinate System is:
PROJCS["NAD83 / UTM zone 13N",

GeoTransform =
319569.9072199312, 4.06320686118055, -0.08151911420909382
4042818.402645736, -0.08151911420909382, -4.06320686118055

Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 319569.907, 4042818.403) (107d 0'53.88"W, 36d30'49.40"N)
Lower Left ( 319215.299, 4025143.453) (107d 0'53.29"W, 36d21'15.87"N)
Upper Right ( 333433.569, 4042540.259) (106d51'36.58"W, 36d30'49.43"N)
Lower Right ( 333078.961, 4024865.310) (106d51'37.13"W, 36d21'15.87"N)
Center ( 326324.434, 4033841.856) (106d56'15.22"W, 36d26' 2.73"N)

Band 1 Block=3412x2 Type=Byte, ColorInterp=Gray

$ gdalwarp test.tif test1.tif
0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo test1.tif
Driver: GTiff/GeoTIFF
Files: test1.tif
Size is 3499, 4418
Coordinate System is:

PROJCS["NAD83 / UTM zone 13N",

Origin = (319215.299073121626861,4042818.402645736001432)
Pixel Size = (4.064024527820449,-4.064024527820449)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 319215.299, 4042818.403) (107d 1' 8.13"W, 36d30'49.16"N)

Lower Left ( 319215.299, 4024863.542) (107d 0'53.06"W, 36d21' 6.79"N)
Upper Right ( 333435.321, 4042818.403) (106d51'36.73"W, 36d30'58.45"N)
Lower Right ( 333435.321, 4024863.542) (106d51'22.84"W, 36d21'16.03"N)
Center ( 326325.310, 4033840.972) (106d56'15.18"W, 36d26' 2.71"N)
Band 1 Block=3499x2 Type=Byte, ColorInterp=Gray



In a GDAL geotransform, the second and last values are not "the pixel size". They are parameters of an affine transformation. However they will be equal to the pixel size where there's no rotation.


The Wikipedia entry on ESRI world files has a good explanation (note GDAL geotransforms are in CABFDE order):




The generic meaning of the six parameters in a world file (as defined by Esri) are:



  • Line 1: A: pixel size in the x-direction in map units/pixel

  • Line 2: D: rotation about y-axis

  • Line 3: B: rotation about x-axis

  • Line 4: E: pixel size in the y-direction in map units, almost always negative

  • Line 5: C: x-coordinate of the center of the upper left pixel

  • Line 6: F: y-coordinate of the center of the upper left pixel


This description is however misleading in that the D and B parameters are not angular rotations, and that the A and E parameters do not correspond to the pixel size if D or B are not zero. The A, D, B and E parameters are sometimes named "x-scale", "y-skew", "x-skew" and "y-scale".



A better description of the A, D, B and E parameters are:



  • Line 1: A: x-component of the pixel width (x-scale)

  • Line 2: D: y-component of the pixel width (y-skew)

  • Line 3: B: x-component of the pixel height (x-skew)

  • Line 4: E: y-component of the pixel height (y-scale), typically negative


All four parameters are expressed in the map units, which are described by the spatial reference system for the raster.


When D or B are non-zero the pixel width is given by:


enter image description here



and the pixel height by:


enter image description here



enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...