在R中相交轮廓和多边形

Posted

技术标签:

【中文标题】在R中相交轮廓和多边形【英文标题】:Intersect the contour and polygon in R 【发布时间】:2013-01-18 05:54:19 【问题描述】:

我已经尝试了几天来创建轮廓,然后在同一个文件上绘制 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 中去除彩色渐变和轮廓,只留下相交部分。

非常感谢任何帮助。

研究:我在 Stack exchange Gis 上找到了这个链接 https://gis.stackexchange.com/questions/25112/clip-depth-contour-with-spatial-polygon,我尝试按照这种方法创建轮廓时总是出错。

我在 https://stat.ethz.ch/pipermail/r-sig-geo/2009-May/005793.html 上发现了另一个类似的帖子。但我无法让它在我的数据集上工作。

我要感谢框中的 Marc 帮助我达到了这一点。

谢谢。

【问题讨论】:

Paul Murrell 在最新的 R-journal 中有一个例子(图 14):journal.r-project.org/archive/2012-2/… 我仍然找不到这个。谁能帮忙???///// 【参考方案1】:

确实,@baptiste 为您提供了解决方案的强烈提示,recent paper by Paul Murrell。 Paul 很慷慨地为我们提供了他整个论文手稿的代码,您可以从his personal website 获得。在副题上,保罗用这篇论文展示了可重复研究的美丽例子。一般来说,我采取了以下方法:

从 shapefile 中提取纬度和经度坐标(执行此操作的函数是 here,由 Paul Hiemstra 编写), 使用您的代码绘制所有内容,

并使用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")

【讨论】:

@Geek On Acid,非常感谢您提供此解决方案。我会查看 Paul 的代码,以便了解更多信息。再次感谢。

以上是关于在R中相交轮廓和多边形的主要内容,如果未能解决你的问题,请参考以下文章

有一个更好的方法吗?

R中多边形和线的复杂裁剪(空间交点?)

增强几何返回相交和相交的不一致结果

多边形相交的简单算法

不相交多边形中的点

如何测试线是否与凸多边形相交?