How to populate gdalwarp with a dynamic value in QGIS geoprocessing?

Question

I am attempting to call gdalwarp (python binding gdal.Warp()) as a script inside QGIS modeler (QGIS 2.18.10, gdal 2.2.0). The goal is to correlate any loaded raster with a UTM zone grid, get the correct EPSG code from that grid, and pipe that code as a dynamic value into gdal.Warp().

To do this, I have written a geoprocessing script because the existing Warp tool does not allow one to dynamically set the EPSG code.

I have looked through the QGIS geoprocessing documentation, but it is a bit thin on examples, which makes troubleshooting difficult.

Here is what I have so far:

##Raster_Layer=raster # Raster to be reprojected
##Joined_Layer=vector # Layer containing the correct EPSG code under the 'EPSG' field
##Output_Raster=output raster

from qgis.core import *
from osgeo import gdal

def get_epsg(lyr):
    for field in lyr.fields():
        fname = field.name()
        if fname == 'EPSG':
            epsg_code = field.value()
    return epsg_code

layer = processing.getObject(Joined_Layer)
srs = get_epsg(layer)

rlayer = processing.getObject(Raster_Layer)
input_raster = gdal.Open(rlayer)
output_raster = Output_Raster
srs_string = 'EPSG:{}'.format(str(srs))
gdal.Warp(output_raster, input_raster, dstSRS=srs_string)

I am sure I'm missing some QGIS-specific wrappers for the api calls, but have no clue where to look for a thorough documentation on using gdal in QGIS geoprocessing scripts (if such a thing exists).

I have a hard time finding good documentation for those calls, but google almost always comes through for me. I googled "python gdal warp set srs" and this was the first hit: https://gis.stackexchange.com/questions/139906/replicating-result-of-gdalwarp-using-gdal-python-bindings You can use that as a template. Note that there are two methods provided: using the API, or, using the command line tool. I rarely use the API, but instead wrap the command line tool in a subprocess() call within my Python scripts. It's easy to set SRS using the command line call: http://www.gdal.org/gdalwarp.html — Jon, Oct 23 '17 at 15:02
Thanks - and I do know how to use the API or CLI in a standalone python script - but what I am looking for are the QGIS wrappers, to make this a user geoprocessing script that can be called into the QGIS modeler. — auslander, Oct 23 '17 at 15:26
Ah, I don't know anything about the QGIS wrappers. I don't see any difference between what you posted and a standalone Python script (i.e. nothing that is QGIS-specific), so probably I don't understand what is wrong with your code. — Jon, Oct 23 '17 at 15:49
Gotcha. things like processing.getObject() and the ##Output_Raster = output raster are particular to QGIS geoprocessing and allow you to write scripts that will offer users a gui tool that accepts inputs or lets them direct output, inside of QGIS. — auslander, Oct 23 '17 at 15:51
When I pull it into the modeler and run it against a raster, it just says "failed: see log"....it reports this message in the log. Which seems redundant and a bit uninformative. Hence the original question :) — auslander, Oct 23 '17 at 15:55
Ah, I see. Does it write "print()" statements to the log? You could make sure the EPSG is being correctly obtained if so. I posted a gdal-based EPSG extractor here: https://gis.stackexchange.com/a/257917/78446 if that is the problem. Otherwise, I'll leave now and let the experts take over. :P — Jon, Oct 23 '17 at 16:05

score 2 · Accepted Answer · answered Oct 24 '17 at 09:53

There's a few minor issues in the script:

epsg_code = field.value() - Assuming each layer contains only one EPSG code in the field, we could use the uniqueValues() method:
```
idx = lyr.fieldNameIndex(fname)
epsg_code = lyr.uniqueValues(idx)[0]
```
gdal.Open() - This requires the path to the raster so we don't need to get its corresponding object. We can replace rlayer = processing.getObject(Raster_Layer) with:
```
rlayer = Raster_Layer
```
You also don't need from qgis.core import *

Here is your slightly modified script which worked for me:

##Raster_Layer=raster
##Joined_Layer=vector
##Output_Raster=output raster

from osgeo import gdal

def get_epsg(lyr):
    for field in lyr.fields():
        fname = field.name()
        if fname == 'EPSG':
            idx = lyr.fieldNameIndex(fname)
            epsg_code = lyr.uniqueValues(idx)[0]
    return epsg_code

layer = processing.getObject(Joined_Layer)
srs = get_epsg(layer)

rlayer = Raster_Layer
input_raster = gdal.Open(rlayer)
output_raster = Output_Raster
srs_string = 'EPSG:{}'.format(str(srs))
gdal.Warp(output_raster, input_raster, dstSRS=srs_string)

Second time in a week you've answered one of my questions, I'd give you double points if I could. This is more than helpful. The overall projection model works end to end and takes less than 1 second to execute. — auslander, Oct 24 '17 at 16:30
@auslander - Glad it helped and thanks for posting this question. I haven't used gdal.warp in this fashion before (always through the processing framework) so I learnt a new method! Also, I don't care for points. Knowledge exchange is the most important reward :) — Joseph, Oct 25 '17 at 09:07

How to populate gdalwarp with a dynamic value in QGIS geoprocessing?

1 Answers1