Как буферизовать векторный шейп-файл с помощью ogr python?

10

Я пытаюсь узнать, как использовать OGR в Python, используя наборы данных о стране и населенных пунктах из http://www.naturalearthdata.com/downloads/50m-cultural-vectors/, Я пытаюсь использовать фильтры и буферы для поиска точек (ne_50m_populated_places.shp) в указанном буфере названной страны (отфильтровано из класса объектов ADMIN в ne_50m_admin_0_countries.shp). Проблема заключается в том, что я не понимаю, какие модули использовать для buffer (). В сценарии я просто использовал произвольное значение 10, чтобы проверить, работает ли сценарий. Скрипт запускается, но возвращает населенные пункты со всего Карибского региона для названной страны «Ангола». В идеале я хочу иметь возможность указать буферное расстояние, скажем, 500 км, но не могу понять, как это сделать, так как я понимаю, что buffer () использует единицы стран.shp, которые будут в формате wgs84 lat / long , Консультации о методе для достижения этой цели будет высоко ценится.

# 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()
user1765941
источник

Ответы:

2

Я могу придумать два варианта:

  1. Вычислите степень, эквивалентную 500 километрам. Затем вы можете ввести его функцию Buffer (). Вы должны быть осторожны, хотя у степени нет постоянного метрического эквивалента. Это зависит от того, на какой широте вы находитесь. Вы можете проверить формулу Haversine, если вы хотите пойти по этому пути.

  2. Другой вариант - проецировать шейп-файл на UTM. Таким образом, вы можете просто использовать 500 километров напрямую. Вы должны найти UTM-зону для вашей области интересов, хотя. (Какой должна быть UTM зона 32S для Анголы, если я не ошибаюсь)

RK
источник
3
Не существует «эквивалента в градусах» ни в 500 километрах, ни на каких-либо других расстояниях, за исключением приблизительно небольших расстояний вблизи экватора: это связано с тем, что соотношение между расстояниями и градусами изменяется как с учетом, так и с широтой. Таким образом, первый вариант обычно не будет работать правильно.
whuber
0
  1. Если вы хотите создать буфер, используя градусы, то учтите, что градусы имеют разное расстояние в зависимости от направления, когда вы не находитесь вблизи экватора. Широта остается неизменной, но долгота на один градус намного меньше в высоких широтах. Ниже моя таблица 500-километровых квадратов в градусах в разных широтах. Я предполагаю, что для Анголы значение 4,4 может быть хорошим предположением, если вам не нужна высокая точность.
  2. Вы можете перепроектировать объекты в python ogr (для этого есть функция Transform) во время чтения, тогда нет необходимости преобразовывать шейп-файл.
500 км при лат. 0,0 составляет 4,491576420597608 x 4,486983030705042 град.
500 км при широте 10,0 составляет 4,491576420597608 x 4,389054945583991 град
500 км при широте 20,0 составляет 4,491576420597608 x 4,16093408959923 град
500 км при широте 30,0 составляет 4,491576420597608 x 3,8117296267699388 град
500 км при широте 40,0 составляет 4,491576420597608 x 3,3535548944407267 град
500 км при широте 50,0 составляет 4,491576420597608 x 2,8010165014556634 град.
500 км при широте 60,0 составляет 4,491576420597608 x 2,170722673038327 град
500 км при широте 70,0 составляет 4,491576420597608 x 1,4808232946314916 град
500 км при широте 80,0 составляет 4,491576420597608 x 0,7505852760718597 град
500 км при широте 84,0 составляет 4,491576420597608 x 0,4516575041056399 град
JaakL
источник
4
Пользователь @Dave X отмечает, что эта таблица ошибочна: фиксированное расстояние охватывает большее количество градусов в более высоких широтах, а не меньше. Похоже, что он мог быть построен путем умножения там, где требуется деление. Тем не менее, это не полностью объясняет расхождения: остаются ошибки порядка нескольких процентов. Как именно вы вычислили эти числа?
whuber