0

I am doing raster by mask on multiple polygons. Selection on the basis of "CONTCD='AFR'" is running fine. Second selection and looping to for mask is only showing all.

layer=QgsVectorLayer("C:/output/b1.shp","b1", "ogr")  
raster = "C:/output/raster.tif" 
output="C:/output/
layer.selectByExpression("CONTCD='AFR'")
countryname=[]
selection = layer.selectedFeatures()
for feat in selection:
    # print an attribute
    countryname.append(feat['ROMNAM'])
    
layer.removeSelection()

for i in countryname:
    # print an attribute
   layer.selectByExpression("ROMNAM"  +'='+ i)
   
   selection = layer.selectedFeatures()
  
   result= gdal.Warp(output + i+".tiff",raster,cutlineDSName=layer)
   layer.removeSelection()`
1
  • Your issue is in gdal.Warp(output + i+".tiff",raster,cutlineDSName=layer). The part cutlineDSName=layer will never work: layer is a QgsvectorLayer and gdal.Warp does not know about this type of object. You are already lucky it does not throw an error.
    – ThomasG77
    Jun 9, 2021 at 9:58

1 Answer 1

0

No need to go through QGIS selection. You can go the "SQL way". There is not dependencies to QGIS: you can run it in the PyQGIS console but throught also Python script outside QGIS if GDAL Python installed

from osgeo import gdal, ogr

output_directory = '/tmp/'
input_raster = 'OB_50M/OB_50M.tif'
input_vector = 'ne_10m_admin_0_countries.shp'
ds_vector = gdal.OpenEx(input_vector)
# Get a subset, similar to your first selection
result_vector_query = ds_vector.ExecuteSQL("SELECT * FROM ne_10m_admin_0_countries WHERE CONTINENT = 'Africa'")

file_prefix = 'whatever_'
# Loop and with cutlineSQL keep the WHERE condition and add your second condition
# Similar to your second selection (in fact, a selection based on a subset)
for i, feature in enumerate(result_vector_query):
    escaped_name = feature.GetField('NAME').replace("'", "\\'")
    print(escaped_name)
    result= gdal.Warp(
        output_directory + file_prefix + str(i) + ".tiff",
        input_raster,
        cutlineDSName=input_vector,
        cutlineSQL=f"SELECT * FROM ne_10m_admin_0_countries WHERE CONTINENT = 'Africa' AND NAME = '{escaped_name}'"
    )

PS: Datasets can be download on NaturalEarthData website to reproduce

1
  • Thank you @ThomasG77.That what i actually want.
    – user94672
    Jun 12, 2021 at 7:35

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.