R код для вырезки куска из shp файла и сохранение его в новом файле
Эта операция особенно актуальна, если приходится обращаться к одному и тому же большому шэйп-файлу много раз, хотя нужен для дела только небольшой его кусок.
##########################
### Cut a part from a large shp file, http://gis.stackexchange.com/questions/54103/how-to-cut-a-shp-based-on-another-shp-using-r
### making polygon that marks the selected area
library(sp) # or other libraries containing sp as a part of them (rgdal, for example)
r4 = cbind(c(6,6,7,7,6), # longitude in the first column
c(54.3,53.7,53.7,54.3,54.3)) # last point in the list of coordinates should be the same as the first one to close the circle
# it has to be column bind (cbind), not the data.frame
r4
sr4=Polygons(list(Polygon(r4)),"r4")
sr_eee=SpatialPolygons(list(sr4), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
sr_eee
plot(depth_20)
plot(sr_eee, add=TRUE)
### end of making polygon for the selected area
#over(sr_eee, depth_20) ### NOT WORKING
### Making intersect between two shapefiles
library(rgeos)
inter<-gIntersection(depth_20, sr_eee) # intersects two shapefiles
inter
str(inter)
plot(inter, col=rgb(0,0,105, maxColorValue=255), axes=TRUE, add=TRUE, border=NA) #plots the result of intersection, but this result is in the format of SpatialPolygons
# to save it in the other shapefile it has to be converted to SpatialPolygonsDataFrame
spp <-SpatialPolygonsDataFrame(inter,data=as.data.frame("inter"))#, proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
spp
plot(spp)
writeOGR(spp, ".", "inter_depth20", driver="ESRI Shapefile") # writes selected intersection between shapefiles into the other shp file, library rgdal is needed
selected_depth20<-readOGR(dsn="C:/_work/Ansgar/ramming_BW_II/dataset from 20130215", layer="inter_depth20") #reads shp file
plot(selected_depth20)
### End of making intersect between two shapefiles
##########################
### Cut a part from a large shp file, http://gis.stackexchange.com/questions/54103/how-to-cut-a-shp-based-on-another-shp-using-r
### making polygon that marks the selected area
library(sp) # or other libraries containing sp as a part of them (rgdal, for example)
r4 = cbind(c(6,6,7,7,6), # longitude in the first column
c(54.3,53.7,53.7,54.3,54.3)) # last point in the list of coordinates should be the same as the first one to close the circle
# it has to be column bind (cbind), not the data.frame
r4
sr4=Polygons(list(Polygon(r4)),"r4")
sr_eee=SpatialPolygons(list(sr4), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
sr_eee
plot(depth_20)
plot(sr_eee, add=TRUE)
### end of making polygon for the selected area
#over(sr_eee, depth_20) ### NOT WORKING
### Making intersect between two shapefiles
library(rgeos)
inter<-gIntersection(depth_20, sr_eee) # intersects two shapefiles
inter
str(inter)
plot(inter, col=rgb(0,0,105, maxColorValue=255), axes=TRUE, add=TRUE, border=NA) #plots the result of intersection, but this result is in the format of SpatialPolygons
# to save it in the other shapefile it has to be converted to SpatialPolygonsDataFrame
spp <-SpatialPolygonsDataFrame(inter,data=as.data.frame("inter"))#, proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
spp
plot(spp)
writeOGR(spp, ".", "inter_depth20", driver="ESRI Shapefile") # writes selected intersection between shapefiles into the other shp file, library rgdal is needed
selected_depth20<-readOGR(dsn="C:/_work/Ansgar/ramming_BW_II/dataset from 20130215", layer="inter_depth20") #reads shp file
plot(selected_depth20)
### End of making intersect between two shapefiles