2

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.

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

1 Answer 1

3

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

3

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.