5

I have a geotiff, and I want to find its extents and reproject them. I can use rasterio:

# Destination CRS
dest_crs = some CRS

# Get projected boundaries with rasterio
clipbox = rasterio.open('mygeotif.tif')
cl, cb, cr, ct = clipbox.bounds
l, b, r, t = transform_bounds(clipbox.crs, dest_crs, clipbox.bounds[0], 
clipbox.bounds[1], clipbox.bounds[2], clipbox.bounds[3])
print(l,b,r,t)

-75.0277857229692 -9.469504586908355 -74.09777685862646 -7.796173999420965

...or I can use GDAL:

# Get projected boundaries with GDAL
clipbox = gdal.Open('mygeotif.tif')
wktproj = clipbox.GetProjection()
sref = osr.SpatialReference()
sref.ImportFromWkt(wktproj)

# Get its extents
ulx, xres, xskew, uly, yskew, yres  = clipbox.GetGeoTransform()
lrx = ulx + (clipbox.RasterXSize * xres)
lry = uly + (clipbox.RasterYSize * yres)

# Reproject the extents
transform = osr.CoordinateTransformation(dest_crf, sref)
tl, tt, garbage = transform.TransformPoint(ulx, uly)
tr, tb, garbage = transform.TransformPoint(lrx, lry)
print(tl, tb, tr, tt)

-75.02766365052823 -9.468344577976415 -74.09777685862646 -7.797125769469675

So one latitude value is the same, but the other three are different. Does anyone have any ideas why this would be the case? I should note that cl == ulx, cb == lry, etc. I.e. the input boundary coordinates are the same pre-projection.

EDIT:

Pre-projection coordinates:

print(ulx,uly,lrx,lry) # GDAL

496950.0 -861870.0 599040.0 -1046760.0

print(cl, ct, cr, cb) # rasterio

496950.0 -861870.0 599040.0 -1046760.0

Projections:

print(dest_crs) # GDAL

"osgeo.osr.SpatialReference; proxy of Swig Object of type 'OSRSpatialReferenceShadow *' at 0x117b85e10"

and using a function called wkt2epsg, we see

print(wkt2epsg(wktproj))

EPSG:32618

print(dest_crs) # rasterio

EPSG:32618

  • 1
    Can you include some of the other values like the bounds before being reprojected, the input and output cell sizes, input/output CRS? – mkennedy Jun 19 '17 at 19:27
  • Added everything except the cell sizes. I'm not sure what you mean since I am only reprojecting four coordinates, not an entire raster. – Jon Jun 19 '17 at 19:50
5

The source coordinate reference system (CRS) is WGS84 / UTM zone 18 North. The points are south of the equator. I'm assuming the destination CRS is WGS84, EPSG:4326.

One of the programs is "reboxing" the raster while the other one isn't. When I unproject all 4 corners of the raster:

496950.0   -861870.0
496950.0  -1046760.0
599040.0  -1046760.0
599040.0   -861870.0

Here's what I get:

-75.02766365052825  -7.797125769469677
-75.02778572296923  -9.469503785651575
-74.09777685862647  -9.468344577976417
-74.10173999154692  -7.796173999420967

Because the raster isn't centered across the central meridian, it's tipped slightly to the northwest. You do not have a perfect rectangle.

rasterio is taking the minimum and maximum values for all 4 corners to return you the actual extent of the output raster. GDAL is returning just the unprojected lower left and upper right coordinates.

  • I don't quite follow. I thought it was: input coordinates -> transforming function -> output coordinates. Are you saying that rasterio or GDAL is doing an extra step? – Jon Jun 19 '17 at 20:30
  • 2
    It certainly looks that way to me. Reboxing is a common problem when reprojecting data. If only two points (of an envelope) are converted, you can easily miss data. The problem could occur with queries too, not just reprojecting a raster or a vector dataset itself. – mkennedy Jun 19 '17 at 20:33
  • Could it have something to do with the transformation assuming coordinates are pixel corners rather than pixel centers? I ask because when I continue the analysis and reproject the full raster using the GDAL output coordinates, the raster is 1/2 pixel off (N-S), but when I use the rasterio coordinates it aligns as expected. – Jon Jun 19 '17 at 20:52
  • 3
    It looks to me that the transform bounds is projecting the envelope and reboxing.. if you project an envelope and return an envelope (pretty much regardless of software) it will be reboxed; the envelope when projected is an irregular quadrilateral, but an envelope must be flat rectangular, so the returned envelope is the min/max of all the X and Y coords of the quadrilateral which will change the width and height slightly. When projecting a point it is unconstrained and will give exactly the location it ends up which is why the transformed box and 2 points give different locations. – Michael Stimson Jun 19 '17 at 21:39
  • I believe you are correct. The source code for transform_bounds shows that it "densifies" the boundaries to account for "non linear transformations." So instead of connecting the corner points with straight lines, it samples reprojected points between them and chooses the max/min as the new envelope. And I think that's what mkennedy was saying to begin with, I just didn't understand. Thanks to all! – Jon Jun 19 '17 at 23:38

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.