12

IPCCプロットで行うように、「重要」な場所を示すために「点描」を追加したいデータがあります

http://www.ipcc.ch/graphics/ar4-wg1/jpg/fig-10-18.jpg

現時点では、R でこれを行うのに本当に苦労しています。

テストデータを作成してプロットすると、次のようになります。

data <- array(runif(12*6), dim=c(12,6) )
over <- ifelse(data > 0.5, 1, 0 )
image(1:12, 1:6, data)

最後にやりたいことは、現在の画像の上にある配列「オーバー」に基づいて、いくつかのポイントをオーバープロットすることです。

助言がありますか！？？

4

4 に答える 4

9

これは役立つはずです - 私は以前に同様のことをして、ここに投稿した関数を書きました。

#required function from www.menugget.blogspot.com
matrix.poly <- function(x, y, z=mat, n=NULL){
 if(missing(z)) stop("Must define matrix 'z'")
 if(missing(n)) stop("Must define at least 1 grid location 'n'")
 if(missing(x)) x <- seq(0,1,,dim(z)[1])
 if(missing(y)) y <- seq(0,1,,dim(z)[2])
 poly <- vector(mode="list", length(n))
 for(i in seq(length(n))){
  ROW <- ((n[i]-1) %% dim(z)[1]) +1
  COL <- ((n[i]-1) %/% dim(z)[1]) +1
 
  dist.left <- (x[ROW]-x[ROW-1])/2
  dist.right <- (x[ROW+1]-x[ROW])/2
  if(ROW==1) dist.left <- dist.right
  if(ROW==dim(z)[1]) dist.right <- dist.left
 
  dist.down <- (y[COL]-y[COL-1])/2
  dist.up <- (y[COL+1]-y[COL])/2
  if(COL==1) dist.down <- dist.up
  if(COL==dim(z)[2]) dist.up <- dist.down
 
  xs <- c(x[ROW]-dist.left, x[ROW]-dist.left, x[ROW]+dist.right, x[ROW]+dist.right)
  ys <- c(y[COL]-dist.down, y[COL]+dist.up, y[COL]+dist.up, y[COL]-dist.down)
  poly[[i]] <- data.frame(x=xs, y=ys)
 }
 return(poly)
}

#make vector of grids for hatching
incl <- which(over==1)

#make polygons for each grid for hatching
polys <- matrix.poly(1:12, 1:6, z=over, n=incl)

    #plot
png("hatched_image.png")
image(1:12, 1:6, data)
for(i in seq(polys)){
    polygon(polys[[i]], density=10, angle=45, border=NA)
    polygon(polys[[i]], density=10, angle=-45, border=NA)
}
box()
dev.off()

ここに画像の説明を入力

または、「点描」の代替：

png("hatched_image2.png")
image(1:12, 1:6, data)
for(i in seq(polys)){
    xran <- range(polys[[i]]$x)
    yran <- range(polys[[i]]$y)
    xs <- seq(xran[1], xran[2],,5)
    ys <- seq(yran[1], yran[2],,5)
    grd <- expand.grid(xs,ys)
    points(grd, pch=19, cex=0.5)
}
box()
dev.off()

ここに画像の説明を入力

アップデート：

Paul Hiemstra のコメントへの (非常に遅い) 応答として、より高い解像度のマトリックスを使用した 2 つの例を次に示します。ハッチングはきれいな規則的なパターンを維持していますが、バラバラになると見栄えがよくありません。点描の例は、はるかに優れています。

n <- 100
x <- 1:n
y <- 1:n
M <- list(x=x, y=y, z=outer(x, y, FUN = function(x,y){x^2 * y * rlnorm(n^2,0,0.2)}))
image(M)
range(M$z)
incl <- which(M$z>5e5)

polys <- matrix.poly(M$x, M$y, z=M$z, n=incl)

png("hatched_image.png", height=5, width=5, units="in", res=400)
op <- par(mar=c(3,3,1,1))
image(M)
for(i in seq(polys)){
  polygon(polys[[i]], density=10, angle=45, border=NA, lwd=0.5)
  polygon(polys[[i]], density=10, angle=-45, border=NA, lwd=0.5)
}
box()
par(op)
dev.off()

ここに画像の説明を入力

png("stippled_image.png", height=5, width=5, units="in", res=400)
op <- par(mar=c(3,3,1,1))
image(M)
grd <- expand.grid(x=x, y=y)
points(grd$x[incl], grd$y[incl], pch=".", cex=1.5)
box()
par(op)
dev.off()

ここに画像の説明を入力

于 2012-07-31T10:19:26.560 に答える
3

これは、ggplot2 を使用した @mdsummer のコメントの精神に基づくソリューションです。最初にグリッドを描画し、次に+特定の値を超えた場所に 'es を描画します。多次元配列や行列ではなく、でggplot2動作することに注意してください。パッケージから使用して、配列/ marix から data.frame フラット構造に変換data.frameできます。meltreshape

geom_tileドキュメントのサンプルデータを使用した具体的な例を次に示します。

pp <- function (n,r=4) { 
 x <- seq(-r*pi, r*pi, len=n) 
 df <- expand.grid(x=x, y=x) 
 df$r <- sqrt(df$x^2 + df$y^2) 
 df$z <- cos(df$r^2)*exp(-df$r/6) 
 df 
} 

require(ggplot2)
dat = pp(200)
over = dat[,c("x","y")]
over$value = with(dat, ifelse(z > 0.5, 1, 0))
ggplot(aes(x = x, y = y), data = dat) + 
   geom_raster(aes(fill = z)) + 
   scale_fill_gradient2() +
   geom_point(data = subset(over, value == 1), shape = "+", size = 1)

ここに画像の説明を入力

于 2012-07-31T10:31:27.980 に答える
3

?image[1]の座標配置メカニズムを使用してそれを行います。

data(volcano)
m <- volcano
dimx <- nrow(m)
dimy <- ncol(m)

d1 <- list(x = seq(0, 1, length = dimx), y = seq(0, 1, length = dimy), z = m)

そのように構築された「画像」を使用すると、オブジェクトの構造とその座標をそのまま維持できます。複数の行列を 3D 配列または複数の要素として収集できますがimage()、それを処理するには拡張する必要があるため、ここではそれらを分けておきます。

データのコピーを作成して、関心のある領域を指定します。

d2 <- d1
d2$z <- d2$z > 155

座標を使用して、関心のあるセルを指定します。非常に大きなラスターがある場合、これはコストがかかりますが、実行するのは非常に簡単です。

pts <- expand.grid(x = d2$x, y = d2$y)
pts$over <- as.vector(d2$z)

プロットを設定します。

op <- par(mfcol = c(2, 1))
image(d1)

image(d1)
points(pts$x[pts$over], pts$y[pts$over], cex = 0.7)

par(op)

ポイントのプロットを変更してさまざまな効果を得ることを忘れないでください。特に、多くのポイントを含む非常に密集したグリッドでは、これらの小さな円をすべて描画するのに時間がかかります。pch = "."良い選択です。

さて、その素敵な投影法にプロットする実際のデータはありますか? 一部のオプションについては、こちらの例を参照してください: http://spatial-analyst.net/wiki/index.php?title=Global_datasets

[1] R には、ラスター データをより高度に処理するためのクラスがあります。2 つの異なるアプローチについては、パッケージ sp および raster を参照してください。

于 2012-07-31T21:45:08.490 に答える
2

これはおそらく遅すぎますが、私の回答も参考として投稿したいと思います。

空間データの便利なオプションの 1 つは、rasterVis パッケージを使用することです。「ベース」ラスター オブジェクトと、点描を描画するために使用する「マスク」オブジェクトを取得したら、次のようなことができます。

require(raster)
require(rasterVis)

# Scratch raster objects
data(volcano)
r1 <- raster(volcano)

# Here we are selecting only values from 160 to 180.
# This will be our "mask" layer.
over <- ifelse(volcano >=160 & volcano <=180, 1, NA) 
r2 <- raster(over)

# And this is the key step:
# Converting the "mask" raster to spatial points
r.mask <- rasterToPoints(r2, spatial=TRUE)

# Plot
levelplot(r1, margin=F) +
layer(sp.points(r.mask, pch=20, cex=0.3, alpha=0.8))

OPが探していたマップに似ています。色、サイズ、タイプなどのポイントのパラメータを微調整できます。?sp.points は、そのために使用できるすべての引数を提供します。

于 2016-05-10T19:55:12.033 に答える