Система координат координат растров и полигонов в 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

Ваш полигональный шейп-файл имеет систему координат, которая не широта - числа очень большие и, возможно, в некоторых системах метры. Назначение строки proj4 не приведет к перепроецированию данных в lat-long, оно просто устанавливает метку того, какие координаты он считает. В этом случае это неправильно!

Вы должны убедиться, что ваши полигоны получают правильную строку proj4 для чисел, которые в них есть - может быть файл [shapefile] .prj вместе с .shp и .dbf, который сообщает вам. Установите для этого строку proj4.

Затем вы можете использовать spTransform из sp или rgdal для проецирования ваших полигонов в координаты WGS84 по длине в широту.

Всегда лучше преобразовывать полигоны в растровые координаты, так как путаница с растровыми координатами может означать перепроецирование всей сетки, которая является грязной ...