Jak buforować wektorowy plik kształtu za pomocą Pythona?

Próbuję nauczyć się korzystać z ogrodu w pythonie przy użyciu zestawu danych o kraju i zaludnionych miejscach z http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. Próbuję użyć filtrów i buforów, aby znaleźć punkty (ne_50m_populated_places.shp) w ​​określonym buforze nazwanego kraju (odfiltrowane z klasy obiektów ADMIN w ne_50m_admin_0_countries.shp). Problem wydaje się polegać na tym, że nie rozumiem, jakich jednostek użyć do buforowania (). W skrypcie po prostu użyłem dowolnej wartości 10, aby sprawdzić, czy skrypt działa. Skrypt działa, ale zwraca zaludnione miejsca z całego regionu Karaibów dla nazwanego kraju „Angola”. Idealnie, chciałbym być w stanie określić odległość bufora, powiedzmy 500 km, ale nie mogę wymyślić, jak to zrobić, ponieważ rozumiem, że bufor () używa jednostek country.shp, który będzie w formacie wgs84 lat / long . Docenione zostaną porady dotyczące sposobu osiągnięcia tego celu.

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

Dwie opcje, o których mogę pomyśleć, to:

1. Oblicz stopień równoważny 500 kilometrów. Następnie możesz wprowadzić funkcję Buffer (). Trzeba być jednak ostrożnym, ponieważ stopień nie ma stałego ekwiwalentu metrycznego. To zależy od szerokości geograficznej. Jeśli chcesz iść tą trasą, możesz sprawdzić formułę Haversine .

2. Inną opcją byłoby rzutowanie pliku kształtu do UTM. W ten sposób możesz bezpośrednio przejechać 500 kilometrów. W strefie zainteresowań znajdziesz strefę UTM. (Która powinna być strefa UTM 32S dla Angoli, jeśli się nie mylę)

1. Jeśli chcesz utworzyć bufor za pomocą stopni, to weź pod uwagę, że stopnie mają całkiem inną odległość w zależności od kierunku, gdy nie jesteś w pobliżu równika. Szerokość pozostaje taka sama, ale jeden stopień długości geograficznej jest znacznie mniejszy na dużych szerokościach geograficznych. Poniżej znajduje się moja tabela z 500-kilometrowymi kwadratami w stopniach na różnych szerokościach geograficznych. Myślę, że dla Angoli wartość 4,4 może być dobra, jeśli nie potrzebujesz wysokiej precyzji.
2. Możesz ponownie rzutować obiekty w Pythonie (jest dla niego funkcja Transform) podczas czytania, wtedy nie ma potrzeby przekształcania plików kształtów.

500 km przy szerokości 0,0 wynosi 4,491576420597608 x 4,486983030705042 st
500 km przy 10.0 lat wynosi 4,491576420597608 x 4,389054945583991 deg
500 km przy lat 20.0 wynosi 4,491576420597608 x 4.16093408959923 st
500 km przy lat 30.0 wynosi 4,491576420597608 x 3,8117296267699388 deg
500 km przy lat 40.0 wynosi 4,491576420597608 x 3,3535548944407267 stopni
500 km przy łatie 50,0 wynosi 4,491576420597608 x 2,8010165014556634 deg
500 km przy łatie 60,0 wynosi 4,491576420597608 x 2,170722673038327 st
500 km przy długości 70,0 wynosi 4,491576420597608 x 1,4808232946314916 deg
500 km przy lat 80.0 wynosi 4,491576420597608 x 0,7505852760718597 deg
500 km w łat 84,0 wynosi 4,491576420597608 x 0,4516575041056399 stopni