[GIS] How to buffer a vector shapefile using ogr python

bufferogrpython

I am trying to learn how to use ogr in python using the country and populated places datasets from http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. I am attempting to use filters and buffers to find points (ne_50m_populated_places.shp) within a specified buffer of a named country (filtered from feature class ADMIN in ne_50m_admin_0_countries.shp). The problem appears to be that I don't understand what units to use for buffer(). In the script I have simply used an arbitrary value of 10 to test if the script works. The script runs but returns populated places from around the Carribean region for named country 'Angola'. Ideally, I want to be able to specify a buffer distance, say 500km, but can't work out how to do this as my understanding is buffer() is using the units of countries.shp that will be in wgs84 lat/long format. Advice on method to achieve this would be much appreciated.

# import modules
import ogr, os, sys


## data source
os.chdir('C:/data/naturalearth/50m_cultural')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open ne_50m_admin_0_countries.shp and get the layer
admin = driver.Open('ne_50m_admin_0_countries.shp')
if admin is None:
  print 'Could not open ne_50m_admin_0_countries.shp'
  sys.exit(1)
adminLayer = admin.GetLayer()

# open ne_50m_populated_places.shp and get the layer
pop = driver.Open('ne_50m_populated_places.shp')
if pop is None:
  print 'could not open ne_50m_populated_places.shp'
  sys.exit(1)
popLayer = pop.GetLayer()

# use an attribute filter to restrict ne_50m_admin_0_countries.shp to "Angola"
adminLayer.SetAttributeFilter("ADMIN = ANGOLA")

# get the Angola geometry and buffer it by 10 units
adminFeature = adminLayer.GetFeature(0)
adminGeom = adminFeature.GetGeometryRef()
bufferGeom = adminGeom.Buffer(10)

# use bufferGeom as a spatial filter on ne_50m_populated_places.shp to get all places
# within 10 units of Angola
popLayer.SetSpatialFilter(bufferGeom)

# loop through the remaining features in ne_50m_populated_places.shp and print their
# id values
popFeature = popLayer.GetNextFeature()
while popFeature:
  print popFeature.GetField('NAME')
  popFeature.Destroy()
  popFeature = popLayer.GetNextFeature()

# close the shapefiles
admin.Destroy()
pop.Destroy()

Best Answer

Two options I can think of are:

  1. Compute the degree equivalent of 500 kilometers. You can then input it Buffer() function. You'd have to be careful though as a degree does not have a constant metric equivalent. It depends on which latitude you're on. You can check out the Haversine formula if you want to go this route.

  2. Another option would be to reproject the shapefile to UTM. That way, you can just use 500 kilometers directly. You'll have find the UTM zone for your area of interest though. (Which should be UTM zone 32S for Angola if I'm not mistaken)

Related Question