1

このrasterパッケージを使用して、ASCII ラスターと ESRI シェープファイルの 2 つのデータセットを読み込みました。ラスターの (水温) データを、湖の海岸線であるシェープファイルの全範囲に抽出したいと考えています。

SpatialPolygonsDataFrame関数を使用して読み込むと、ESRI シェープファイルは として扱われshapefile()ます。

shapefile <- shapefile("shore.shp",verbose=TRUE)

関数を使用raster()して ASCII ラスターを読み込みました。

raster <- raster("1995_001.asc",native=TRUE,crs="+proj=merc +datum=WGS84 +ellps=WGS84 +units=m")

シェープファイルの座標参照情報は次のとおりです。

+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0

ラスターのそれ (つまりcrs、関数内の引数を使用して次のように強制されraster()ます):

+proj=merc +datum=WGS84 +ellps=WGS84 +units=m +towgs84=0,0,0

次にspTransform()、パッケージ内の関数を使用しrgdalて、シェープファイルの空間参照をラスターの空間参照に強制しました。

spTransform(shapefile, CRS(projection(raster)))

最後に、以下を提出しました。

extract(raster,shapefile,method="simple",fun=mean,small=TRUE,na.rm=TRUE,df=FALSE)

ただし、型のオブジェクトをextract()返すだけです。この問題は、座標参照の明示的な強制から生じていると思います。NULLlist

さらに、show()各データセットで関数を使用した結果は次のとおりです。

> show(raster) class : RasterLayer dimensions : 1024, 1024, 1048576 (nrow, ncol, ncell) resolution : 1800, 1800 (x, y) extent : -10288022, -8444822, 4675974, 6519174 (xmin, xmax, ymin, ymax) coord. ref. : NA data source : in memory names : layer values : -9999, 8.97 (min, max)

> show(shapefile) class : SpatialPolygonsDataFrame features : 1 extent : 597568.5, 998261.6, 278635.3, 668182.2 (xmin, xmax, ymin, ymax) coord. ref. : +proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 variables : 3 names : AREA, PERIMETER, HECTARES min values : 59682523455.695, 5543510.075, 5968252.346 max values : 59682523455.695, 5543510.075, 5968252.346

これらのフォーラムで解決策のない同様の質問を多数検索しました。誰か（仮想の）手を貸してくれませんか？

事前にどうもありがとうございました。

4

1 に答える 1

1

これ (空間的に重複していないように見えるデータ) は、座標参照系に関する混乱から生じる一般的な問題です。これは、いくつかの理由で発生する可能性があります。

1) 領域は実際には重なっていません

2) RasterLayer への crs の割り当てが間違っている可能性があります (またはポリゴンの割り当て)。オブジェクト ( show(raster)) を表示してください。rasterこれは、パッケージについて質問するときに、ほとんどの場合役に立ちます。

3) spTransform の結果を代入しません。Rは通常、「インプレース」変更を行わないことに注意してください。必ず行うshapefile <- spTransform(shapefile, crs(raster))

常に健全性テストを行い、問題が解決するまで元に戻してください。ここで開始する場所は、次のことです。

plot(raster)
plot(shapefile, add=TRUE)

ラスター データの上にポリゴンがプロットされているかどうかを確認します。

それでもうまくいかない場合は、エクステントを見てください。たとえば、次のようになります (ちなみに、これは、他の人が実際に回答できる、自己完結型の再現可能な例/質問を作成する方法も示しています):

library(raster)
# the spatial parameters of your raster
r <- raster(ncol=100, nrow=100, xmn=-10288022, xmx=-8444822, ymn=4675974, ymx=6519174, crs="+proj=merc +datum=WGS84")
values(r) <- 1:ncell(r)

# the extent of your SpatialPolygons
ep <- extent(597568.5, 998261.6, 278635.3, 668182.2)
p <- as(ep, 'SpatialPolygons')
crs(p) <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83"

# tranform both data set to longlat
rgeo <- projectRaster(r, crs='+proj=longlat +datum=WGS84')
pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

# plot on top of Google map
library(dismo)
e <- union(extent(rgeo), extent(pgeo))
g <- gmap(e,lonlat=TRUE)
plot(g, inter=TRUE)
plot(extent(rgeo), add=TRUE, col='red', lwd=2)
plot(pgeo, add=TRUE, col='blue', lwd=2)

ここに画像の説明を入力

明らかに、2 つのデータ ソースは重複していません。どちらが正しいかはあなただけが知っているので、もう一方が間違った場所にある理由を調べ始めることができます。

于 2015-06-05T00:34:58.400 に答える