[GIS] GeoTIFF projections with Python GDAL

coordinate systemgdalgeoreferencingpyprojpython

I have a GeoTIFF I imported into Python environment with GDAL. The image has spatial reference data, which are the parameters describing the geo model:

(-104.31292762744124, 0.00030593847623094916, 0.0, 32.77290675637778, 0.0, -0.0002583980094641447)

There is no Projection metadata ( raster.GetProjection() ). How can I use these coordinates to create a projection (I want real world coordinates for each pixel location)?

I looked into osr.SpatialReference and pyproj but I couldn't find any simple or quick explanations. Sorry if this is basic but the online documentation just wasn't helpful for me.

Best Answer

That is the affine transformation matrix, which describes (in this order in your case): [x coordinate of top left corner, pixel width, rotation on x axis, y coordinate of top left corner, rotation on y axis, pixel height]. Fortunately, this is enough to do what you want.

You might find it a little easier to use rasterio, which is an interface to GDAL that is much more user friendly than the GDAL Python bindings.

Here is how it would work:

# you import rasterio like this
import rasterio

# you can open a dataset like this
with rasterio.open("myfile.tif") as dataset:

   # you can see the affine transformation matrix (that you had above) like this:
   print(dataset.transform)

   # you can convert from pixel coordinates to geographical like this:
   # (this example uses the middle of the dataset)
   x, y = dataset.xy(dataset.height // 2, dataset.width // 2)

   # and vice-versa like this:
   # (replace x and y with some coordinates that fall within your dataset)
   row, col = dataset.index(x, y)

Note that Image Space (pixel) coordinates are in the order y,x (row, col), whereas coordinate space uses x, y (lng lat / easting northing). It is also worth knowing that Image Space has the origin in the top left, where Coordinate Space has the orgin in the bottom left, (hence the pixel height being negative in your dataset).

If you want to add the CRS as well (not necessary for the above) - you could do this:

   # assign CRS as WGS84
   dataset.crs = rasterio.crs.CRS.from_string('EPSG:4326')

I have guessed your Coordinate Reference System (CRS) there - if it is a different one to WGS84, then you can just get the correct EPSG code from here and swap it out with 4236 in the above snippet.

Alternatively - if you didn't want to switch to use rasterio, you could do convert from pixel locations (px_x, px_y) to geographical coordinates using your affine transformation matrix (t) like this:

# convert px_x and px_y rto geographical coordinates
x, y = (topleft_x + (px_x * resolution_x), topleft_y - (px_y * resolution_x))

# so in your case, if t is the affine transformation object you had:
x, y = (t[0] + (px_x * t[1]), t[3] - (px_y * t[5]))
Related Question