Bagaimana cara buffer shapefile vektor menggunakan ogr python?

10

Saya mencoba mempelajari cara menggunakan ogr di python menggunakan negara dan mengisi kumpulan data dari http://www.naturalearthdata.com/downloads/50m-cultural-vectors/. Saya mencoba menggunakan filter dan buffer untuk menemukan poin (ne_50m_populated_places.shp) dalam buffer yang ditentukan dari negara yang disebutkan (difilter dari ADMIN kelas fitur di ne_50m_admin_0_countries.shp). Masalahnya tampaknya saya tidak mengerti unit apa yang digunakan untuk buffer (). Dalam skrip saya hanya menggunakan nilai arbitrer 10 untuk menguji apakah skrip berfungsi. Script berjalan tetapi mengembalikan tempat-tempat berpenduduk dari seluruh wilayah Karibia untuk negara bernama 'Angola'. Idealnya, saya ingin dapat menentukan jarak buffer, katakanlah 500km, tetapi tidak bisa mengetahui bagaimana melakukan ini karena pemahaman saya adalah buffer () menggunakan unit country.shp yang akan berada dalam format wgs84 lat / panjang . Saran tentang metode untuk mencapai ini akan sangat dihargai.

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

Jawaban:

2

Dua opsi yang dapat saya pikirkan adalah:

  1. Hitung derajat yang setara dengan 500 kilometer. Anda kemudian dapat memasukkannya fungsi Buffer (). Anda harus berhati-hati karena gelar tidak memiliki metrik yang setara. Itu tergantung pada lintang mana Anda berada. Anda dapat melihat formula Haversine jika Anda ingin pergi rute ini.

  2. Pilihan lain adalah untuk r memproyeksikan shapefile ke UTM. Dengan begitu, Anda bisa langsung menggunakan 500 kilometer. Anda akan menemukan zona UTM untuk bidang yang Anda minati. (Yang seharusnya menjadi zona UTM 32S untuk Angola jika saya tidak salah)

RK
sumber
3
Tidak ada "derajat setara" 500 kilometer, atau jarak lain, kecuali sekitar jarak kecil dekat khatulistiwa: itu karena hubungan antara jarak dan derajat berubah dengan bantalan serta lintang. Jadi opsi pertama biasanya tidak akan berfungsi dengan benar.
whuber
0
  1. Jika Anda ingin membuat buffer menggunakan Derajat, maka pertimbangkan bahwa derajat memiliki jarak yang sangat berbeda tergantung arah ketika Anda tidak berada di dekat khatulistiwa. Lintang tetap sama, tetapi satu derajat bujur jauh lebih kecil di lintang tinggi. Di bawah ini adalah meja saya seluas 500 km persegi dalam derajat di lintang berbeda. Saya kira untuk Angola nilai 4,4 mungkin merupakan tebakan yang baik jika Anda tidak membutuhkan presisi tinggi.
  2. Anda dapat memproyeksikan ulang objek dalam python ogr (ada fungsi Transform untuknya) selama membaca, maka tidak perlu membuat konversi konversi.
500 km pada lat 0,0 adalah 4,491576420597608 x 4,486983030705042 deg
500 km pada lat 10.0 adalah 4.491576420597608 x 4.389054945583991 deg
500 km di lat 20.0 adalah 4.491576420597608 x 4.16093408959923 deg
500 km pada lat 30.0 adalah 4.491576420597608 x 3.8117296267699388 deg
500 km di lat 40.0 adalah 4.491576420597608 x 3.3535548944407267 deg
500 km di lat 50.0 adalah 4.491576420597608 x 2.8010165014556634 deg
500 km di lat 60.0 adalah 4.491576420597608 x 2.170722673038327 deg
500 km di lat 70.0 adalah 4.491576420597608 x 1.4808232946314916 deg
500 km di lat 80.0 adalah 4.491576420597608 x 0.7505852760718597 deg
500 km di lat 84.0 adalah 4.491576420597608 x 0.4516575041056399 deg
JaakL
sumber
4
Pengguna @Dave X mencatat bahwa tabel ini salah: jarak tetap mencakup lebih banyak derajat pada garis lintang lebih tinggi, bukan lebih rendah. Tampaknya itu mungkin dibangun dengan membuat perkalian di mana divisi diperlukan. Meski begitu, itu tidak sepenuhnya menjelaskan perbedaan: masih ada kesalahan pada urutan beberapa persen. Bagaimana tepatnya Anda menghitung angka-angka ini?
whuber