Rのラスター&ポリゴンの調整参照システム
-
26-10-2019 - |
質問
私にはいくつかのポリゴンがあり、これらのポリゴン内のいくつかのラスター層から平均値を抽出するのが好きです。それらをARCMAPに追加したとき、2つのデータ型の予測が一致しないことに気付きました。プロジェクトツール(データ管理ツールボックス>プロジェクションと変換ツールセット>ラスター)を使用して、ArcGISのディスプレイの問題を解決できます。そこで、次の方法でデータをRにロードすることにより、プロジェクションを標準化しようとしました(コードの一部)。
ラスター:
for (i in 1:length(rasterlist1))
{ndvi_raster_stack1[i]<-raster(rasterlist1[i])
raster::NAvalue(ndvi_raster_stack1[[i]])<--999
projection(ndvi_raster_stack1[[i]])<-"+proj=utm +ellps=WGS84 +datum=WGS84 +units=m"}
> ndvi_raster_stack1[[1]]
class : RasterLayer
dimensions : 226, 150, 33900 (nrow, ncol, ncell)
resolution : 0.57504, 0.5753628 (x, y)
extent : -28.728, 57.528, -55.08, 74.952 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0
values : Z:\master\lusmeg_sw_kernel_data\ndvi0910\Y2008_P47.tif
min value : -91
max value : 550.8125
ポリゴン:
for (i in 1:length(poplist))
{pop_kernels[i]<-readShapeSpatial(poplist[i],repair=TRUE,proj4string=CRS("+proj=utm +ellps=WGS84 +datum=WGS84 +units=m"))
pop_kernels[[i]]<-unionSpatialPolygons(pop_kernels[[i]],ID=c(rep(1,times=length(pop_kernels[[i]])-1),0),threshold=NULL,avoidGEOS=FALSE)}
> str(pop_kernels[[1]])
Formal class 'SpatialPolygons' [package "sp"] with 4 slots
..@ polygons :List of 2
.. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. .. .. ..@ Polygons :List of 2
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 2404422 893343
.. .. .. .. .. .. ..@ area : num 1.15e+12
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:1625, 1:2] 2551236 2533236 2533236 2523236 2523236 ...
.. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 2468549 865776
.. .. .. .. .. .. ..@ area : num 6.31e+11
.. .. .. .. .. .. ..@ hole : logi TRUE
.. .. .. .. .. .. ..@ ringDir: int -1
.. .. .. .. .. .. ..@ coords : num [1:1385, 1:2] 2551236 2551236 2563236 2563236 2569236 ...
.. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. ..@ plotOrder: int [1:2] 1 2
.. .. .. ..@ labpt : num [1:2] 2404422 893343
.. .. .. ..@ ID : chr "0"
.. .. .. ..@ area : num 1.15e+12
.. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
.. .. .. ..@ Polygons :List of 1
.. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
.. .. .. .. .. .. ..@ labpt : num [1:2] 2468549 865776
.. .. .. .. .. .. ..@ area : num 6.31e+11
.. .. .. .. .. .. ..@ hole : logi FALSE
.. .. .. .. .. .. ..@ ringDir: int 1
.. .. .. .. .. .. ..@ coords : num [1:1385, 1:2] 2551236 2541236 2541236 2529236 2529236 ...
.. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2
.. .. .. .. .. .. .. .. ..$ : NULL
.. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y"
.. .. .. ..@ plotOrder: int 1
.. .. .. ..@ labpt : num [1:2] 2468549 865776
.. .. .. ..@ ID : chr "1"
.. .. .. ..@ area : num 6.31e+11
..@ plotOrder : int [1:2] 1 2
..@ bbox : num [1:2, 1:2] 1819236 207017 3013236 1577017
.. ..- attr(*, "dimnames")=List of 2
.. .. ..$ : chr [1:2] "x" "y"
.. .. ..$ : chr [1:2] "min" "max"
..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
.. .. ..@ projargs: chr " +proj=utm +ellps=WGS84 +datum=WGS84 +units=m +towgs84=0,0,0"
ポリゴンとラスターを個別にプロットできますが、ポリゴンの1つをラスター上にプロットしようとすると、ポリゴンは表示されません。
plot(ndvi_raster_stack1[[1]],xlab="Longitude",ylab="Latitude")
plot(pop_kernels[[1]],col="black",add=TRUE)
彼らはまだ「オーバーラップ」していないようです。これは、さまざまな境界ボックスによっても示されています。
> bbox(ndvi_raster_stack1[[1]])
min max
s1 -28.728 57.528
s2 -55.080 74.952
> bbox(pop_kernels[[1]])
min max
x 1819236 3013236
y 207017 1577017
ポリゴン内のラスター値を抽出したいので、それらが正しい方法で参照されることを確認する必要があります。誰かが問題が何であるかを考えていますか?
解決
Polygon ShapeFileには、LATのない座標系があります。数字は非常に大きく、おそらく一部のシステムではメートルです。 proj4Stringを割り当てることは、データをLAT-LONGに再損傷することはありません。これは、自分が考えている調整のラベルを設定するだけです。この場合、それは間違っています!
あなたのポリゴンが彼らが持っている数値に対して正しいproj4stringを取得することを確認する必要があります - あなたに伝える.shpおよび.dbfとともに[ShapeFile] .Prjファイルがあるかもしれません。 proj4stringをそれに設定します。
次に、SPまたはRGDALのsptransformを使用して、ポリゴンをLAT LONG WGS84座標に投影できます。
ポリゴンをラスター座標に変換するのが常に最善です。ラスター座標をいじることは、グリッド全体を非難することを意味する可能性があるため、これは乱雑です...