R中柵格和多邊形的座標參考系統

Question

我有幾個多邊形&我喜歡從這些多邊形內的幾個柵格圖層中提取平均值。 當我將這些添加到ArcMap時，我意識到兩種數據類型的投影不匹配。我可以通過使用項目工具（數據管理工具箱>投影和轉換工具集>柵格）來解決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" 

我可以繪製多邊形和柵格分開，但是當我嘗試在柵格上繪製其中一個多邊形時，多邊形不會顯示：

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 

因爲我想多邊形中提取柵格值，我必須確保他們以正確的方式引用。 有人有一個想法可能是什麼問題？

Answer 1

您的多邊形shapefile的座標系不經緯 - 數字非常大，可能在某些系統中爲米。分配一個proj4string不會將數據重新映射到lat-long，它只是設置它認爲它是什麼座標的標籤。在這種情況下，它的錯誤！

你需要確保你的多邊形得到正確的proj4string的數字，他們有他們 - 可能有一個[shapefile] .prj文件以及.shp和.dbf告訴你。將proj4string設置爲。

然後，您可以使用sp或rgdal的spTransform將您的多邊形投影到經緯度爲WGS84的座標。

它總是最好的將多邊形轉換爲柵格座標，因爲搞亂柵格座標可能意味着重新投影整個柵格，這是凌亂的...