QGIS Python Console – Extracting Raw Data from WCS Raster Using QGIS Python Console

clippyqgisqgis-python-consolerasterwcs

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.

Best Answer

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

Related Question