Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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()
share|improve this question

2 Answers 2

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)

share|improve this answer
1  
There is no "degree equivalent" of 500 kilometers, nor of any other distance, except approximately for small distances near the equator: that's because the relationship between distances and degrees changes with bearing as well as latitude. Thus the first option will usually not work correctly. –  whuber Feb 28 '13 at 20:41
  1. If you want to make buffer using Degrees, then consider that degrees have quite different distance depending on direction when you are not near equator. Latitude stays same, but one degree longitude is much smaller in high latitudes. Below is my table of 500-km squares in degrees in different latitudes. I guess that for Angola value of 4.4 may be good guess if you do not need high precision.
  2. You can reproject objects in python ogr (there is Transform function for it) during reading, then there is no need to shapefile conversions.
500 km at lat 0.0 is 4.491576420597608 x 4.486983030705042 deg
500 km at lat 10.0 is 4.491576420597608 x 4.389054945583991 deg
500 km at lat 20.0 is 4.491576420597608 x 4.16093408959923 deg
500 km at lat 30.0 is 4.491576420597608 x 3.8117296267699388 deg
500 km at lat 40.0 is 4.491576420597608 x 3.3535548944407267 deg
500 km at lat 50.0 is 4.491576420597608 x 2.8010165014556634 deg
500 km at lat 60.0 is 4.491576420597608 x 2.170722673038327 deg
500 km at lat 70.0 is 4.491576420597608 x 1.4808232946314916 deg
500 km at lat 80.0 is 4.491576420597608 x 0.7505852760718597 deg
500 km at lat 84.0 is 4.491576420597608 x 0.4516575041056399 deg
share|improve this answer
1  
User @Dave X notes that this table is erroneous: a fixed distance spans a greater number of degrees at higher latitudes, not a lesser. It appears it might have been constructed by making a multiplication where a division is required. Even so, that does not fully explain the discrepancies: there remain errors on the order of several percent. Exactly how did you compute these numbers? –  whuber Feb 28 '13 at 20:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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