r - Rで輪郭とポリゴンを交差させます

Question

数日間、輪郭を作成してから、シェープファイルと輪郭を同じファイルにプロットしようとしています。これで、同じプロット上に等高線とシェープファイルを作成できるようになりました。シェープファイルで輪郭をクリップしたいのですが、シェープファイルのみを表示します。

データtemp.csvはこのリンクにありますhttps://www.dropbox.com/s/mg2bo4rcr6n3dks/temp.csvシェープファイル は次の場所にあります：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)

上記のコードの出力は次のとおりです。

ここで、ポリゴンシェープファイルから色付きのグラデーションと輪郭を削除し、交差部分のみを残します。

どんな助けでも大歓迎です。

調査：スタック交換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で別の同様のスレッドを見つけました。しかし、データセットで機能させることができませんでした。

この点に到達するのを手伝ってくれたボックスのマークに感謝します。

ありがとう。

Answer 1

確かに、@ baptisteは、解決策の強力なヒントを提供しました。これは、PaulMurrellによる最近の論文です。ポールは、彼の個人的なウェブサイトから入手できる彼の紙の原稿全体のコードを私たちに提供してくれました。サイドトピックでは、ポールはこの論文で再現可能な研究の美しい例を示しています。一般的に、私は次のアプローチを取りました。

• シェープファイルから緯度と経度の座標を抽出します（これを行う関数は、Paul Hiemstraによるものです）。
• コードですべてをプロットし、

• polypath抽出された座標をベースラインとして使用して、シェープファイルで定義された境界の外側にあるすべてのものを削除するために使用します。

  #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")