4

I am working on a script that will take an unprojected raster, and using a UTM grid (shapefile), will obtain the correct SRS (I have a shapefile where each zone is attributed with its numerical SRS) and perform a warp/projection using GDAL on that raster. It is a quicker way than figuring out the right zone, and then manually doing the warp.

Procedurally, I am grabbing the raster corners and constructing WKT geometry (polygon); the objective is to do a spatial join with the correct feature in the UTM grid. Here's what I have so far:

import os
import ogr
from osgeo import gdal

def getRasterBox(raster):
    src = gdal.Open(raster)
    ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform()
    lrx = ulx + (src.RasterXSize * xres)
    lry = uly + (src.RasterYSize * yres)
    return (uls, uly, lrx, lry)

def findUTMzone(raster_box, utmgrid):
    # do the stuff

utmzone_path = r'C:\temp'
utmzone_file = 'utm_grid.shp'
utmzone = os.path.join(utmzone_path, utmzone_file)

utm_ds = ogr.Open(utmzone)
utm_layer = utm_ds.GetLayer()

rasterdir = r'C:\temp'
raster_file = 'n38.tif'
raster = os.path.join(rasterdir, raster_file)

rasterdetails = getRasterBox(raster)
rasterWKT = 'POLYGON(({} {}, {} {}, {} {}, {} {}, {} {}))'.format(rasterdetails[0], rasterdetails[1], rasterdetails[0], rasterdetails[3], rasterdetails[2], rasterdetails[3], rasterdetails[2], rasterdetails[1], rasterdetails[0], rasterdetails[1])

Using the QuickWKT plugin for QGIS, I can see that I'm getting a good WKT ring for the raster bounds. However, I'm unsure how to perform the spatial join in a way that grabs the feature attribute containing the correct SRS code. How do I do this part (under def findUTMzone(raster_box, utmgrid):)? I have seen a few arcpy examples, but am looking for something open-source.

  • Are you open to a QGIS python solution? – artwork21 Oct 18 '17 at 17:53
  • Not married to one approach. I can easily change the interpreter for the project to pyqgis, but would need a little extra help on the glue code, as i'm not familiar with the qgis.core package. – auslander Oct 18 '17 at 18:01
  • Using Q the general logic would be to search for "pyqgis standalone app", then "pyqgis create layer or in-memory layer from bbox coordinates", and finally "pyqgis spatial join". If you pull all this logic together and update your question and let us know where you are stuck I'm sure someone will be able to help you out. – artwork21 Oct 18 '17 at 18:46
  • Perhaps the better question is: how do I do a spatial join with GDAL/OGR? Is such a thing possible? Currently I have no problems creating a geometry from bbox. I just want to take geometry A and spatial join it to shapefile b. – auslander Oct 18 '17 at 18:49
  • Have you seen this post, spatialthoughts.com/blog/gis/ogr-spatial-join – artwork21 Oct 18 '17 at 18:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.