1

I have a polygon layer and I want to know the altitude of the centroids. I guess I have to find on the internet (how?) a layer that contains data about the territory's height (in this case Italy in particular Toscana and Liguria) and to find a way to correlate it with my layer (using Python is better). Does anyone have a suggestion?

vlayer = QgsProject.instance().mapLayersByName('Centroids')[0]
rlayer = QgsProject.instance().mapLayersByName('Elevation')[0]
params = {'INPUT_VECTOR': QgsProcessingFeatureSourceDefinition(vlayer.id(), True), 
          'INPUT_RASTER': QgsProcessingFeatureSourceDefinition(rlayer.id(), True), 
          'COLUMN_PREFIX' : 'Altitude',
          'OUTPUT': 'memory:'}

newLayer = processing.run("qgis:rastersampling", params)
layerOutput = newLayer['OUTPUT']
QgsProject.instance().addMapLayer(layerOutput)

I tried in this way, but there's something wrong. This is the error message: Could not load source layer for INPUT: no value specified for parameter

Improve this question
2
  • 3
    You will need a Digital Elevation Model (DEM). You can download tiles for your area of interest here (90m): dwtkns.com/srtm or here (30m): dwtkns.com/srtm30m. Create centroids and run 'Sample raster values' algorithm, then join elevation attribute back to polygon layer. If you want a Python solution you should try to make a code attempt.
    – Ben W
    Commented May 15, 2020 at 0:27
  • Thanks both. I did it, now I'm trying to do the same task using python. The only way I could think is using processing.run algorithm ("qgis:rastersampling") but i'm not sure how to set params and I can't find any example.
    – Lorenzo
    Commented May 15, 2020 at 13:29

1 Answer 1

Reset to default
2

Try:

vlayer = QgsProject.instance().mapLayersByName('New scratch layer')[0]
rlayer = QgsProject.instance().mapLayersByName('nh_62_3')[0]

centroids = [g.geometry().centroid() for g in vlayer.getFeatures()]

#Create a point layer from polygon centroids
vl = QgsVectorLayer("Point?crs={}&index=yes".format(vlayer.crs().authid()), "myLayer", "memory")
provider = vl.dataProvider()
for centroid in centroids:
    f = QgsFeature()
    f.setGeometry(centroid)
    provider.addFeature(f)

#Set z from raster
pointz = processing.run("native:setzfromraster", 
    {'INPUT':vl,
    'RASTER':rlayer.source(),'BAND':1,'NODATA':0,'SCALE':1,
    'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']

#Extract z value
pointz2 = processing.run("native:extractzvalues", 
    {'INPUT':pointz,'SUMMARIES':[0],'COLUMN_PREFIX':'z_',
    'OUTPUT':'TEMPORARY_OUTPUT'})['OUTPUT']

QgsProject.instance().addMapLayer(pointz2)

enter image description here

Improve this answer
2
  • Thanks it works. Is there a way to maintain in the new layer pointz2 the attributes of the first vector layer?
    – Lorenzo
    Commented May 16, 2020 at 13:38
  • Spatial join / join attributes by location. Or start with centroid that has attributes.
    – Bera
    Commented May 16, 2020 at 13:41

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.