Membaca dan Menulis Shapefile dengan R

Format Esri Shapefile merupakan salahs atu format yang populer untuk penyimpanan data spasial berbasis vektor seperti garis, titik dan poligon. pada kesempatan kali ini akan diuraikan sedikit mengenai bekerja menggunakan software R untuk membaca dan menulis format Shapefile. Terdapat beberapa packages/plugin di R yang dapat digunakan untuk membaca dan menulis format SHP.

1. Menggunakan rgdal

merupakan salah satu paket R yang berfungsi memberikan interface GDAL/OGR library ke software R, yang memungkinkan import/export data geospasial. plugin tersebut termasuk fungsi readOGR dan writeOGR yang bukan saja bisa membaca dan menulis file SHP tapi juga beberapa jenis file vektor. Dengan tambahan ogrInfo yang berguna untuk mendapatkan detail dari file tanpa harus membaca keseluruhan dataset. Fungsi ini juga bisa secara otomatis membaca dan menulis proyeksi jika ada.

2.Menggunakan Maptools

Paket maptools termasuk di dalamnya beberapa fungsi berguna seperti membaca, menulis, konversi, dan menghandel obyek spasial di R. fungsi umum untuk membaca dan menulis pada shapefile adalah readShapeSpatial dan writeSpatialShape, kedua fungsi ini biasanya secara otomatis langsung mengenali apakah shp atau obyek R mempunyai entitas titik, garis atau poligon, dan kemudian membacanya atau/dan menulisnya menggunakan fungsi yang lebih khusus lagi. Fungsi khusus seperti readShapeLines untuk membaca garis, juga dapat dipanggil secara langsung. Manfaat pemanggilan langsung adalah error akan segera muncul jika anda menggunakan tipe data yang keliru. Tidak seperti gdal, packages maptools tidak membaca dan menulis informasi proyeksi, jadi anda harus mengolahnya secara manual.

3. Menggunakan PBSMapping

Paket ini dapat membaca format SHP tetapi tidak dapat menulisnya. tetapi perlu dicatat bahwa fungsi PBSMapping menggunakan tipe data spasial kustom yang dioptimasi untuk bekerja dengan berbagai macam paket fungsi, sehingga sangat sulit untuk dapat mengambil keuntungan dari fungsi yang didefinisikan oleh paket dari sp, walaupun paket maptools menyediakan fungsi konversi ke dalam format berbeda.

CODE R :

rgdal :

library(rgdal)

# for shapefiles, first argument of the read/write/info functions is the
# directory location, and the second is the file name without suffix

# optionally report shapefile details
ogrInfo(“.”, “nw-rivers”)
# Source: “.”, layer: “nw-rivers”
# Driver: ESRI Shapefile number of rows 12
# Feature type: wkbLineString with 2 dimensions
# +proj=longlat +datum=WGS84 +no_defs
# Number of fields: 2
# name type length typeName
# 1 NAME 4 80 String
# 2 SYSTEM 4 80 String

# read in shapefiles
centroids.rg rivers.rgcounties.rg

# note that readOGR will read the .prj file if it exists
print(proj4string(counties.rg))
# [1] ” +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0″

# generate a simple map showing all three layers
plot(counties.rg, axes=TRUE, border=”gray”)
points(centroids.rg, pch=20, cex=0.8)
lines(rivers.rg, col=”blue”, lwd=2.0)

# write out a new shapefile (including .prj component)
writeOGR(counties.rg, “.”, “counties-rgdal”, driver=”ESRI Shapefile”)

maptools :

library(maptools)

# read in shapefiles; here we use the specialized readShape* functions,
# but readShapeSpatial would produce identical output in all three cases
centroids.mp rivers.mpcounties.mp

# note that readShape* does _not_ read the shapefile’s .prj file
print(proj4string(counties.mp))
## [1] NA

# specifying projection information is not strictly necessary for
# plotting, but does yield nicer default axis labels and aspect ratio in
# the case of geographic data
proj4string(counties.mp) <- “+proj=longlat +datum=WGS84″

# generate a simple map showing all three layers
plot(counties.mp, axes=TRUE, border=”gray”)
points(centroids.mp, pch=20, cex=0.8)
lines(rivers.mp, col=”blue”, lwd=2.0)

# write out a new shapefile (but without .prj); the more general
# writeSpatialShape would produce equivalent output
writePolyShape(counties.mp, “counties-maptools”)

PBSMapping :

library(PBSmapping)

# read in shapefiles
centroids.pb rivers.pbcounties.pb

# note that importShapefile reads the .prj file if it exists, but it
# does not adopt the proj4 format used by the above approaches
proj.abbr proj.full print(proj.abbr)
# [1] “LL”

# generate map using PBSmapping plotting functions
plotPolys(counties.pb, projection=proj.abbr, border=”gray”,
xlab=”Longitude”, ylab=”Latitude”)
addPoints(centroids.pb, pch=20, cex=0.8)
addLines(rivers.pb, col=”blue”, lwd=2.0)

CONTOH DATA DAN KODE DAPAT DIDOWNLOAD DI SINI

sumber : http://www.nceas.ucsb.edu/

Tinggalkan Pesan