7

I am working with a shapefile containing the district polygons of the political map of India and a raster file containing the elevation data of the same country. My goal is to compute the average and variance altitude within 10 km from the borders of each district. Using Zonal Statistics I was able to compute the statistics I am looking for at the district level, the further step I need is disregard the pixels lying in the central part.

I have tried to get rid of the internal parts of the districts with add rings, but the user should draw the internal polygons by himself. This is not what I am looking for: I want to have an internal polygon with exactly the same angles as the original one but with smaller sides, computed automatically for all India.

How can I create, alternatively, a mask in the shapefile for the area that I don't want (or that I want, the one that most fits the case) to take into account and then use it to compute mean and variance with the Zonal Statistics plugin?

The ideal answer would provide me a PyQGIS code but it's fine also using icons in QGIS.

  • Have you tried buffering the districts with a -10 km distance in order to get your new zones? – Tom Mar 7 '16 at 16:24
  • I've removed ArcGIS from your question to enable it to be focused on PyQGIS which is your stated preference. – PolyGeo Mar 7 '16 at 20:37
3

I tried out the next code with a polygon shapefile (as active layer) where, first, it is produced a "negative" buffer (-10000 m) as memory layer. Afterward, this memory layer is used as difference layer for the active layer for producing a polygon without the internal part (it is obtained a central ring).

layer = iface.activeLayer()

feat = layer.getFeatures().next()

#determining buffer -10000 m
buffer = feat.geometry().buffer(-10000,-1)

#Extract CRS from layer
CRS = layer.crs().postgisSrid()

URI = "Polygon?crs=epsg:"+str(CRS)+"&field=id:integer""&index=yes"

#Create polygon layer for buffer
mem_layer = QgsVectorLayer(URI,
                           "buffer",
                           "memory")

#Prepare mem_layer for editing
mem_layer.startEditing()

#Set feature for buffer
feat2 = QgsFeature()

#Set geometry for buffer
feat2.setGeometry(buffer)

#set attributes values for buffer
feat2.setAttributes([1])

mem_layer.addFeature(feat2, True)

#stop editing and save changes
mem_layer.commitChanges()

#preparing layer for difference
feat3 = mem_layer.getFeatures().next()

difference = feat.geometry().difference(feat3.geometry())

#Create polygon layer for difference
mem_layer2 = QgsVectorLayer(URI,
                           "diference",
                           "memory")

#Prepare mem_layer2 for editing
mem_layer2.startEditing()

#Set feature for difference
feat2 = QgsFeature()

#Set geometry for difference
feat2.setGeometry(difference)

#set attributes values for difference
feat2.setAttributes([1])

mem_layer2.addFeature(feat2, True)

#stop editing and save changes
mem_layer2.commitChanges()

#add Map Layer to Registry
QgsMapLayerRegistry.instance().addMapLayer(mem_layer2)

Next image was obtained after running the code at Python Console of QGIS:

enter image description here

It works as expected.

Editing Note:

Next version of the script works for all features of a polygon vector layer.

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

#determine buffer -10000 m
buffer = [ feat.geometry().buffer(-10000,-1) for feat in feats ] 

#Extract CRS from layer
CRS = layer.crs().postgisSrid()

URI = "Polygon?crs=epsg:"+str(CRS)+"&field=id:integer""&index=yes"

#Create polygon layer for buffer
mem_layer = QgsVectorLayer(URI,
                           "buffer",
                           "memory")

#Prepare mem_layer for editing
mem_layer.startEditing()

n = len(feats)

#Set feature for buffer
feats2 = [ QgsFeature() for i in range(n) ]

for i in range(n): 
    #Set geometry for buffer
    feats2[i].setGeometry(buffer[i])

    #set attributes values for buffer
    feats2[i].setAttributes([i])

    mem_layer.addFeature(feats2[i], True)

#stop editing and save changes
mem_layer.commitChanges()

feats3 = [ feat3 for feat3 in mem_layer.getFeatures() ]

difference = []

for i in range(n):
    difference.append(feats[i].geometry().difference(feats3[i].geometry()))

#Create polygon layer for difference
mem_layer2 = QgsVectorLayer(URI,
                           "difference",
                           "memory")

#Prepare mem_layer for editing
mem_layer2.startEditing()

#Set feature for difference
feats2 = [ QgsFeature() for i in range(n) ] 

for i in range(n):
    #Set geometry for difference
    feats2[i].setGeometry(difference[i])

    #set attributes values for difference
    feats2[i].setAttributes([i])

    mem_layer2.addFeature(feats2[i], True)

#stop editing and save changes
mem_layer2.commitChanges()

#add Map Layer to Registry
QgsMapLayerRegistry.instance().addMapLayer(mem_layer2)

It was run with a vector layer where, in some features, the distance from the borders is lesser than 10 km; as it can be observed at the next image:

enter image description here

  • The figure you posted above is exactly what I wanted, however your script it doesn't run correctly in my case: first it does not produce holes in the various geometries, second it returns only the first district of India living outside the other 592. I have set the projection as India 1975 UTM but nothing change with respect the 4326 WGS 84 projection. The first district is an archipelagos in which not all islands are bigger than 20 km in all directions. Can this be a problem? – Paolo Santini Mar 8 '16 at 14:12
  • I modified the script. It works for all features. The archipelagos will not be a problem. – xunilk Mar 8 '16 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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