1

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.

1

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]))
| improve this answer | |
0

Assuming you know the spatial reference of your data (I'd guess it is WGS84 by taking a look at the origin coordinates you provided and the pixel size), you can set the spatial reference this way:

from osgeo import gdal
from osgeo import osr

# create WGS84 Spatial Reference
sr = osr.SpatialReference()
sr.ImportFromEPSG(4326)

# open raster in writing mode
fn = r'C:\path\to\my_raster.tif'
ds = gdal.Open(fn, 1)

# set spatial reference
ds.SetProjection(sr.ExportToWkt())

# save raster
ds = None
| 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.