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()
