Como armazenar buffer em um shapefile de vetor usando ogr python?

10

Estou tentando aprender como usar ogr em python usando o país e conjuntos de dados de locais populosos em http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. Estou tentando usar filtros e buffers para encontrar pontos (ne_50m_populated_places.shp) dentro de um buffer especificado de um país nomeado (filtrado da classe de recurso ADMIN em ne_50m_admin_0_countries.shp). O problema parece ser que eu não entendo quais unidades usar para buffer (). No script, simplesmente usei um valor arbitrário de 10 para testar se o script funciona. O script é executado, mas retorna lugares populosos da região do Caribe para o país nomeado 'Angola'. Idealmente, eu quero poder especificar uma distância de buffer, digamos 500 km, mas não consigo descobrir como fazer isso, pois meu entendimento é que buffer () está usando as unidades de countries.shp que estarão no formato wgs84 lat / long . Os conselhos sobre o método para conseguir isso seriam muito apreciados.

# 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
fonte

Respostas:

2

Duas opções em que posso pensar são:

  1. Calcule o equivalente em grau de 500 quilômetros. Você pode então inserir a função Buffer (). Você teria que ter cuidado, pois, um diploma não tem um equivalente métrico constante. Depende da latitude em que você está. Você pode verificar a fórmula Haversine se quiser seguir esse caminho.

  2. Outra opção seria projetar novamente o shapefile no UTM. Dessa forma, você pode simplesmente usar 500 quilômetros diretamente. Você encontrará a zona UTM para sua área de interesse. (Qual deve ser a zona UTM 32S para Angola, se não me engano)

RK
fonte
3
Não há "equivalente em graus" de 500 quilômetros, nem em qualquer outra distância, exceto aproximadamente em pequenas distâncias próximas ao equador: isso ocorre porque a relação entre distâncias e graus muda com o rumo e com a latitude. Portanto, a primeira opção geralmente não funcionará corretamente.
whuber
0
  1. Se você deseja fazer um buffer usando Graus, considere que os graus têm uma distância bastante diferente, dependendo da direção, quando você não está perto do equador. A latitude permanece a mesma, mas a longitude de um grau é muito menor em altas latitudes. Abaixo está minha tabela de quadrados de 500 km em graus em diferentes latitudes. Acho que para Angola o valor de 4,4 pode ser um bom palpite se você não precisar de alta precisão.
  2. Você pode reprojetar objetos no python ogr (existe a função Transform) durante a leitura e, portanto, não há necessidade de moldar as conversões.
500 km a lat 0,0 0.0 é 4.491576420597608 x 4.486983030705042 deg
500 km a lat 10.0 é 4.491576420597608 x 4.389054945583991 deg
500 km a 20.0 lat é 4.491576420597608 x 4.16093408959923 deg
500 km a lat 30.0 é 4.491576420597608 x 3.8117296267699388 deg
500 km a lat. 40.0 é 4.491576420597608 x 3.3535548944407267 deg
500 km a lat 50.0 é 4.491576420597608 x 2.8010165014556634 deg
500 km a lat 60.0 é 4.491576420597608 x 2.170722673038327 deg
500 km a lat. 70.0 é 4.491576420597608 x 1.4808232946314916 deg
500 km a lat 80.0 é 4.491576420597608 x 0.7505852760718597 deg
500 km a lat 84.0 é 4.491576420597608 x 0.4516575041056399 deg
JaakL
fonte
4
O usuário @Dave X observa que esta tabela está incorreta: uma distância fixa abrange um número maior de graus em latitudes mais altas, e não menor. Parece que pode ter sido construído fazendo uma multiplicação onde uma divisão é necessária. Mesmo assim, isso não explica totalmente as discrepâncias: permanecem erros na ordem de vários por cento. Exatamente como você calculou esses números?
whuber