Memilih nilai yang benar untuk proj4string untuk pembacaan shapefile di R?

14

Saya memiliki shapefile poligon dan file CSV lain yang berisi daftar titik sebagai (Lat, Lng) berpasangan ..

Saya ingin memeriksa setiap pasangan (lat, lng) dari file CSV yang poligonnya berada di dalam ..

Shapefile diproyeksikan dan file proj berbunyi seperti ini:

PROJCS["Transverse_Mercator",GEOGCS["GCS_OSGB 1936",
DATUM["D_OSGB_1936",SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",49],PARAMETER["central_meridian",-2],PARAMETER["scale_factor",0.9996012717],PARAMETER["false_easting",400000],PARAMETER["false_northing",-100000],UNIT["Meter",1]]

Rencana saya adalah sebagai berikut:

  1. Baca shapefile menggunakan readShapePolyfungsi dalam MapToolspaket R.
  2. Baca titik koordinat dari file CSV ke dalam dataframe dan konversikan ke SpatialPointsDataFrame
  3. Gunakan overfungsi untuk menentukan poligon yang berada di dalamnya.

Untuk melakukannya, saya perlu menentukan proj4stringsaat memuat shapefile di langkah 1 dan juga mengubah koordinat dari file CSV ke sistem proyeksi yang sama menggunakan spTransformfungsi sebelum menerapkan overfungsi di langkah 3 karena mengharuskan titik dan poligon harus berada di bawah sistem proyeksi yang sama.

Adakah ide tentang nilai apa yang seharusnya untuk konten file proj yang ditunjukkan di atas?

coordinate-system shapefile r proj Moustafa Alzantot
sumber
Jika shapefile Anda sudah memiliki proyeksi yang ditentukan, gunakan "readOGR" dalam paket rgdal. Paket ini adalah pembungkus untuk GDAL dan benar-benar menggantikan fungsi baca / tulis shapefile di maptools. Fungsi ini menangani semua jenis topologi dan mempertahankan informasi proyeksi.
Jeffrey Evans
Ketika saya mencoba memuat file bentuk menggunakan readOGRfungsi saya selalu mendapatkan Tidak dapat membuka kesalahan file
Moustafa Alzantot
OK, Sekarang saya sudah bisa membaca file menggunakan readOGR .. menggunakan summaryfungsi untuk SpatialPolygonDataFrameobjek memberi saya nilai yang benar untukproj4string
Moustafa Alzantot
Nah, tanpa perincian tentang bagaimana Anda menggunakan fungsi itu sulit untuk membantu Anda! Bagian dari sintaks adalah direktori tempat data berada dan Anda tidak memerlukan ekstensi .shp dalam nama file. Sesuatu seperti readOGR (getwd (), "YourShape") harus berfungsi jika Anda memiliki direktori-kerja Anda disetel ke tempat yang sama dengan shepfile Anda.
Jeffrey Evans
Terima kasih @JeffreyEvans, itu berhasil sekarang dan saya menggunakannya untuk mendapatkan proj4string
Moustafa Alzantot

Jawaban:

14

Proj4string adalah string PROJ4 crs yang valid .

lihat Bagaimana saya bisa mendapatkan string proj4 atau kode EPSG dari file .prj shapefile? dan Shapefile PRJ ke PostGIS SRID lookup table?

pendeknya:

  • Anda dapat menggunakan gdalinfo seperti pada referensi pertama atau binding Python GDAL seperti pada referensi kedua.

Atau

  • buka Prj2EPSG (layanan sederhana untuk mengonversi informasi proyeksi teks yang terkenal dari file .prj menjadi kode EPSG standar)
  • Masukkan konten file prj Anda

masukkan deskripsi gambar di sini

  • hasilnya adalah EPSG: 27700 jadi versi pertama dari string PROJ4 adalah

    " + init = epsg: 27700 "

`Atau

masukkan deskripsi gambar di sini

  • klik pada Proj4 dan string PROJ4 yang lengkap adalah:

    " + proj = tmerc + lat_0 = 49 + lon_0 = -2 + k = 0,9996012717 + x_0 = 400000 + y_0 = -100000 + ellps = lapang + datum = OSGB36 + unit = m + no_defs "

gen
sumber
10

Berikut adalah situs web yang sangat berguna untuk mengambil kode EPSG untuk proyeksi yang diberikan. Dalam kasus Anda, proyeksi adalah "EPSG: 27700". Jika Anda memiliki proyeksi yang ditentukan untuk shapefile Anda dapat menetapkan proyeksi ketika Anda membuat SpatialPointsDataFrame dan kemudian menggunakan definisi proyeksi dari shapefile yang Anda impor. Menggunakan "readOGR" dari paket rgdal akan mempertahankan definisi proyeksi. Berikut ini adalah contoh penetapan dan penarikan string koordinat pada data kelas sp.

require(sp)
require(rgdal)

# Use meuse dataset
data(meuse)

# Coerce into SpatialPointsDataframe
coordinates(meuse) <- ~x+y

# Assign projection
proj4string(meuse) <- CRS("+init=epsg:28992")

# Pull some observations and transform to Lat/Long
meuse.geo <- meuse[sample(dim(meuse)[1],10),]
  prj.LatLong <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
    meuse.geo <- spTransform(meuse.geo, prj.LatLong)

# Pull projection string from meuse.geo and use in spTransform
#   to reproject meuse to lat/long  
( prj <- proj4string(meuse.geo) )   
meuse <- spTransform(meuse, CRS(prj))
Jeffrey Evans
sumber