I am trying to learn how to use ogr in python using the country and populated places datasets from http://www.naturalearthdata.com/downloads/10m-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 script runs but returns populated places from around the Carribean region. 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. 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 a better 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)
    angolaGeom = adminFeature.GetGeometryRef()
    bufferGeom = angolaGeom.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()