Sistema di coordinate di riferimento di Raster e poligoni in R

Question

Ho diversi poligoni e mi piace per estrarre i valori medi da diversi strati raster all'interno di questi poligoni. Quando ho aggiunto quelli di ArcMap mi sono reso conto che le proiezioni dei due tipi di dati non corrispondono. Ho potuto risolvere il problema per la visualizzazione in ArcGIS utilizzando lo strumento del progetto (toolbox Data Management> Proiezioni e trasformazioni set di strumenti> Raster). Così ho cercato di standardizzare la proiezione caricando i dati nel R nel modo seguente (parte del codice):

raster:

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

I poligoni:

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"

posso tracciare i poligoni e le raster separatamente, ma quando provo a tracciare uno dei poligoni oltre un raster, i poligoni non vengono visualizzati:

plot(ndvi_raster_stack1[[1]],xlab="Longitude",ylab="Latitude")
plot(pop_kernels[[1]],col="black",add=TRUE)

Sembra che ancora non si "sovrappongono". Questo è anche indicato dai diversi riquadri di limitazione, penso:

> 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

Perché voglio estrarre i valori raster all'interno dei poligoni, devo essere sicuro che essi si fa riferimento in modo corretto. Qualcuno ha un'idea di che cosa il problema potrebbe essere?

Answer 1

Your polygon shapefile has a coordinate system that isnt lat-long - the numbers are very large and probably metres in some system. Assigning a proj4string won't reproject the data to lat-long, it just sets the label of what coordinates it thinks it is. In this case, its wrong!

You need to make sure your polygons get the right proj4string for the numbers they have in them - there may be a [shapefile].prj file along with the .shp and .dbf that tells you. Set the proj4string to that.

Then you can use spTransform from sp or rgdal to project your polygons to lat-long WGS84 coordinates.

Its always best to transform polygons to raster coordinates, since messing with raster coordinates can mean reprojecting the whole grid, which is messy...