3

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).

Improve this question
7
  • 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: gis.stackexchange.com/questions/139906/… 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: gdal.org/gdalwarp.html Commented Oct 23, 2017 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. Commented Oct 23, 2017 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. Commented Oct 23, 2017 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. Commented Oct 23, 2017 at 15:51
  • Pardon my denseness, but why doesn't your script work? Commented Oct 23, 2017 at 15:52

1 Answer 1

Reset to default
2

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)
Improve this answer
2
  • 1
    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. Commented Oct 24, 2017 at 16:30
  • 1
    @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 :) Commented Oct 25, 2017 at 9:07

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.