1

I have a DEM file and shapefile of specific area, and I want to Create a module in Python that writes name of each area which is shown in the picture, with the greatest differences in altitude on the territory (min-max range of the altitude).

I also have to create a new layer containing the shape of these areas, together with all its original attributes.

I have GDAL and ArcPy libraries.

Could you please help me determine what I should do?

closed as off-topic by PolyGeo Dec 10 '15 at 1:30

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "Questions seeking help to debug/write/improve code must include the desired behavior, a specific problem or error and the shortest code necessary to reproduce it in the question itself. Providing a clear problem statement and evidence of a code attempt will help others to help you. See: How to create a Minimal, Complete, and Verifiable example." – PolyGeo
If this question can be reworded to fit the rules in the help center, please edit the question.

  • What GIS library for Python are you using? – PolyGeo Dec 4 '15 at 23:44
  • i am using gdal library. – Saeed Dec 5 '15 at 0:25
  • There's some sample code on the GDAL/OGR cookbook. If I understand you right, I think you're asking about 'zonal statistics'. – Steven Kay Dec 5 '15 at 0:30
  • yes this is about zonal statistics. – Saeed Dec 5 '15 at 0:47
  • Have you tried the sample code suggested by @StevenKay? If so, and if you are still stuck, then include a code snippet to show what you tried. – PolyGeo Dec 6 '15 at 2:48
2

If it doesn't have to be Python, just use R, it's very well suited for this with high-level tools.

Pseudo-session (requires rgdal installed as well for read/write shapefile):

library(raster)
s <- shapefile("file.shp")
d <- raster("dem.whatever")
s$demmin <- extract(d, s, fun = min, na.rm = TRUE)
s$demmax <- extract(d, s, fun = max, na.rm = TRUE)

Write it back out

shapefile(s, "file_minmax.shp")
1

As it was pointed out by Steven Kay, there's some sample code on the GDAL/OGR cookbook about 'zonal statistics', so it's some extensive when it is compared with the mdsumner's R code. I know that R power can be used in Python by using the rpy and rpy2 python modules. After install them, I tried out next code in PyQGIS environment:

from PyQt4.QtCore import *

import rpy2.robjects as robjects
import os

r = robjects.r

mapcanvas = iface.mapCanvas()
layers = mapcanvas.layers()

vp = layers[0].dataProvider()
vpath_tmp = vp.dataSourceUri().split('|')
vpath = vpath_tmp[0]

rp = layers[1].dataProvider()
rpath = rp.dataSourceUri()

vroot, vfile = os.path.split(vpath)
rroot, rfile = os.path.split(rpath)

r.setwd(vroot)

r.library('raster')

arg1 = 's <- shapefile(' + "\"" + vfile + "\"" + ')'
arg2 = 'd <- raster(' + "\"" + rfile + "\"" + ')'

r(arg1)
r(arg2)

print "Wait..."

demmin = r('s$demmin <- extract(d, s, fun = min, na.rm = TRUE)')
demmax = r('s$demmax <- extract(d, s, fun = max, na.rm = TRUE)')

fields = [ QgsField('min', QVariant.Double), QgsField('max', QVariant.Double) ]
vp.addAttributes( fields )
layers[0].updateFields()

idx1 = layers[0].fieldNameIndex('min')
idx2 = layers[0].fieldNameIndex('max')

n = len(demmin)

for i in range(n):
    new_values = { idx1 : demmin[i], idx2 : demmax[i] }
    vp.changeAttributeValues( {i:new_values} )

print "Done!"

with the polygon vector layer and dem raster of next image:

enter image description here

After running the code at the Python Console of QGIS, two new fields (min and max) were incorporated to attributes table of vector layer and preserved its original attributes; as it can be observed at next image:

enter image description here

It worked for me.

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