我已经尝试了几天来创建轮廓,然后在同一个文件上绘制 shapefile 和轮廓。现在,我可以在同一个图上创建轮廓和 shapefile。我想用 shapefile 裁剪轮廓,只显示 shapefile。
数据 temp.csv 可以在此链接上找到 https://www.dropbox.com/s/mg2bo4rcr6n3dks/temp.csv
Shapefile 可以在以下位置找到:https://www.dropbox.com/sh/ztvmibsslr9ocmc/YOtiwB8p9p
脚本文件 image.scale.R 可以在以下位置找到“https://www.dropbox.com/s/2f5s7cc02fpozk7/image.scale.R”
到目前为止,我使用的代码如下:
## Required packages
library(maptools)
library(rgdal)
library(sp)
library(maptools)
library(sm)
require(akima)
require(spplot)
library(raster)
library(rgeos)
## Set Working Directory
setwd("C:\\Users\\jdbaba\\Documents\\R working folder\\shape")
## Read Data from a file
age2100 <- read.table("temp.csv",header=TRUE,sep=",")
x <- age2100$x
y <- age2100$y
z <- age2100$z
####################################
##Load the shape file
#####################################
shapefile <- readShapePoly("Export_Output_4.shp")
fld <- interp(x,y,z)
par(mar=c(5,5,1,1)) filled.contour(fld)
###Import the image.scale
source source("image.scale.R")
# http://menugget.blogspot.de/2011/08/adding-scale-to-image-plot.html
x11(width=8, height=7)
layout(matrix(c(1,2), nrow=1, ncol=2), widths=c(6,1), height=6, respect=TRUE)
layout.show(2)
par(mar=c(4,4,1,2))
image(fld,axes=T)
contour(fld, add=TRUE)
#points(age2100$x,age2100$y, pch=".", cex=2,legend=F)
plot(shapefile,add=T,lwd=2)
box()
par(mar=c(4,0,1,4))
image.scale(fld, xlab="Eastings", ylab="Northings", xaxt="n", yaxt="n", horiz=FALSE)
axis(4)
mtext("Salinity", side=4, line=2.5)
上述代码的输出如下:
现在,我想从多边形 shapefile 中去除彩色渐变和轮廓,只留下相交部分。
任何帮助都受到高度赞赏。
研究:我找到了这个链接 https://gis.stackexchange.com/questions/25112/clip-depth-contour-with-spatial-polygon在堆栈交换 Gis 上,我尝试遵循此方法,但在创建轮廓时总是出错。
我在 https://stat.ethz.ch/pipermail/r-sig-geo/2009-May/005793.html 上发现了另一个类似的帖子.但我无法让它在我的数据集上工作。
我想在框中感谢 Marc 帮助我达到这一点。
谢谢。
最佳答案
事实上,@baptiste 为您提供了解决方案的强烈提示,recent paper by Paul Murrell .保罗慷慨地向我们提供了他整个论文手稿的代码,您可以从 his personal website 获得代码。 .在附带的主题上,保罗用这篇论文展示了可重复研究的漂亮例子。一般来说,我采取了以下方法:
polypath
使用提取的坐标作为基线,删除 shapefile 定义的边界之外的所有内容。#function to extract coordinates from shapefile (by Paul Hiemstra)
allcoordinates_lapply = function(x) {
polys = x@polygons
return(do.call("rbind", lapply(polys, function(pp) {
do.call("rbind", lapply(pp@Polygons, coordinates))
})))
}
q = allcoordinates_lapply(shapefile)
#extract subset of coordinates, otherwise strange line connections occur...
lat = q[110:600,1]
long = q[110:600,2]
#define ranges for polypath
xrange <- range(lat, na.rm=TRUE)
yrange <- range(long, na.rm=TRUE)
xbox <- xrange + c(-20000, 20000)
ybox <- yrange + c(-20000, 20000)
#plot your stuff
plot(shapefile, lwd=2)
image(fld, axes=F, add=T)
contour(fld, add=T)
#and here is the magic
polypath(c(lat, NA, c(xbox, rev(xbox))),
c(long, NA, rep(ybox, each=2)),
col="white", rule="evenodd")
关于r - 在 R 中与轮廓和多边形相交,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/14393172/