1

I am trying to replicate the following manual workflow in the Python console of QGIS 3.4.7

Manual workflow:

1) Connect to WCS service.

Name: USGS National Map Elevation OGC-WCS

URL: https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer?request=GetCapabilities&service=WCS

2) Add the elevation data to the map canvas. I did this by dragging in the WCS layer "DEP3ElevationPrototype"

3) Zoomed to my area of interest

4) Right click the WCS layer > Export > Save As

5) Checked the box for Raw Data. Specified output file name and extent. Clicked OK.

Problem
The above workflow worked perfectly for me. However, I would like to have this same workflow in a python script. I couldn't find a processing tool for "Export > Save As", so I tried using GDAL Clip raster by mask layer

import processing
from qgis.core import *

### Raster layer from WCS
rlayername = 'elev'
uri = QgsDataSourceUri()
uri.setParam('url', 'https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer?request=GetCapabilities&service=WCS')
uri.setParam("identifier", rlayername)
rlayer = QgsRasterLayer(str(uri.encodedUri()), 'DEP3ElevationPrototype', 'wcs')


mask_filepath = 'N:\zTEMP\Ian\QGIS_Scripting\Black Diamon Study Area.shp'
mask_layer = QgsVectorLayer(mask_filepath, 'mask', 'ogr')


### Params for GDAL Clip Raster by Mask
parameters = {'INPUT': rlayer,
            'MASK': mask_layer,
            'NODATA': -9999,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': 'N:\zTEMP\Ian\QGIS_Scripting\DEP3_Extract_Clip.tif'}


### Clip
print ("running algorithm")
processing.runAndLoadResults("gdal:cliprasterbymasklayer",parameters)
print("done")

I assume the problem is that Clip Raster by Mask cannot take a WCS raster, though I could be mistaken. Are there any ways to doing something like the manual workflow in the Python console?

Note that in this case I don't really care about the distinction clipping by extent vs clipping by mask. In the end I just want to get a TIF extract from a WCS raster.

  • Can't you use WCS subsetting to get the aoi? – nmtoken Apr 27 at 12:26
  • you have an error in the wcs itself, that ideintifier does not exist. See the capablities – Fran Raga Apr 29 at 8:23
  • My QGIS 2.18.13 says "AttributeError: 'QgsRasterProjector' object has no attribute 'setCrs'" – Hugo Jul 4 at 13:44
1

this algorithm doesn't accept WCS so you have to look for an alternative.

Note: I correct the load of the WCS, since its code doesn't work for me

I propose this solution

import processing
from qgis.core import *
import urllib.parse

# Make WCS Uri
def makeWCSuri( url, layer ):  
    params = {  'dpiMode': 7 ,
                'identifier': layer,
                'url': url.split('?')[0]  } 

    uri = urllib.parse.unquote( urllib.parse.urlencode(params)  )
    return uri 

### Raster layer from WCS
rlayername = 'DEP3ElevationPrototype'
wcsUri = makeWCSuri('https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer', rlayername )
rlayer = QgsRasterLayer(wcsUri, 'DEP3ElevationPrototype', 'wcs')
if not rlayer.isValid():
  print ("Layer failed to load!")

#QgsProject.instance().addMapLayer(rlayer)

mask_filepath = r'D:\temp\mask2.shp'
mask_layer = QgsVectorLayer(mask_filepath, 'mask', 'ogr')
#QgsProject.instance().addMapLayer(mask_layer)

# Save raster
renderer = rlayer.renderer()
provider = rlayer.dataProvider()
crs = rlayer.crs()

pipe = QgsRasterPipe()
projector = QgsRasterProjector()
projector.setCrs(provider.crs(), provider.crs())

if not pipe.set(provider.clone()):
    print("Cannot set pipe provider")

# Commented for extract raw data
# if not pipe.set(renderer.clone()):
    # print("Cannot set pipe renderer")

if not pipe.insert(2, projector):
    print("Cannot set pipe projector")

out_file = 'D:/temp/temporal.tif'
file_writer = QgsRasterFileWriter(out_file)
file_writer.Mode(1)

print ("Saving")

extent = mask_layer.extent()

opts = ["COMPRESS=LZW"]
file_writer.setCreateOptions(opts)
error = file_writer.writeRaster(
    pipe,
    extent.width (),
    extent.height(),
    extent,
    crs)

if error == QgsRasterFileWriter.NoError:
    print ("Raster was saved successfully!")
    layer = QgsRasterLayer(out_file, "result")
    QgsProject.instance().addMapLayer(layer)
else:
    print ("Raster was not saved!")

I hope it helps

  • I've run your script and I was able to get a geotiff with the correct extent. However, the output raster appears to contain the RGB values rather than the raw elevation data (i.e. the raster consists only of pixels RGB = 17,17,17). I consulted the documentation for the QGsRasterFileWriter class and changed the Mode to 0 (raw). However, the output was the same as when Mode =1. Is there anything else I need to modify to get the raw values? – isevans77 Apr 29 at 15:07
  • updated the answer to get the raw information – Fran Raga Apr 29 at 16:52

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.