Como rasterizar SpatialPolygons em R?

10

Estou tentando extrair os valores de batimetria da minha área de interesse de uma camada raster de batimetria mundial usando a função 'rasterize' no pacote {sp}.

* Edições: Encontrei a função 'extrair', que parece ser mais o que estou procurando.

Isto é o que eu fiz até agora:

> class(subarea0) #This is my area of interest (Eastern Canadian Arctic Sea)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"

> extent(subarea0)
class       : Extent 
xmin        : -82.21997 
xmax        : -57.21667 
ymin        : 60.2 
ymax        : 78.16666

library(marelac)
data("Bathymetry")#World bathymetric data in library (marelac)
names(Bathymetry);class(Bathymetry);str(Bathymetry)
[1] "x" "y" "z"
[1] "list"
List of 3
 $ x: num [1:359] -180 -179 -178 -177 -176 ...
 $ y: num [1:180] -89.5 -88.5 -87.5 -86.5 -85.5 ...
 $ z: num [1:359, 1:180] 2853 2873 2873 2873 2873 ...

  raster_bath<-raster(Bathymetry)#Transformed into a raster layer
    extent(raster_bath) <- extent(subarea0)#Transform the extend of my raster to the extend of my SpatialPolygons

>ras_sub0<-rasterize(subarea0,raster_bath)#rasterize my SpatialPolygons (*Edits: not the function that I need here, but I am still interested to learn what results mean)
Found 2 region(s) and 10 polygon(s)
> plot(ras_sub0)
> plot(subarea0, add=TRUE)

insira a descrição da imagem aqui

> ras_sub0
class       : RasterLayer 
dimensions  : 180, 359, 64620  (nrow, ncol, ncell)
resolution  : 0.06964709, 0.0998148  (x, y)
extent      : -82.21997, -57.21667, 60.2, 78.16666  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 
values      : in memory
min value   : 1 
max value   : 2 
layer name  : layer 

Eu não entendo o resultado. Por que estou recebendo duas cores para cada um dos meus polígonos? O que eles querem dizer?

Como posso finalmente obter o contorno da profundidade batimétrica? Isso tem algo a ver com minha resolução ou alteração das dimensões?

* Edições: ok, agora fiz o seguinte:

v <- extract(raster_bath, subarea0)#Extract data from my Raster Layer for the locations of my SpatialPolygons

v é uma lista e não tenho certeza de como / sob que forma religar essas informações com meu polígono espacial ...

Obrigado!

GodinA
fonte
Quando você diz 'obter contorno de profundidade de batimetria', quer dizer que também deseja gerar contornos?
Simbamangu

Respostas:

6

Sua linha ras_sub0<-rasterize(subarea0,raster_bath)está apenas pegando o número do índice dos polígonos e atribuindo-o aos valores da varredura.

Se você deseja apenas a interseção do seu polígono e da varredura:

subarea0_bathy <- intersect(raster_bath, subarea0)

Atualização : como o @GodinA observa, parece que intersect () às vezes não retorna uma varredura que tem toda a extensão do polígono! Para contornar isso, você pode cruzar com uma varredura um pouco maior que a original:

r2 <- raster() # create a generic raster
extent(r2) <- extent(subarea0) + 1 # add 1 degree of extent (0.5 each side)
r3 <- intersect(raster_bath, r2) 
Simbamangu
fonte
Obrigado @Simbamangu, tentei o comando "interseção", no entanto, quando plogo a subárea0_bathy e meu polígono espacial (subárea0), eles não se sobrepõem completamente, por que? Para responder à sua pergunta anterior, sim, eu gostaria de eventualmente também ter meu contorno batimétrico. Alguma sugestão?
Godina
Sim, olhei para a saída dos dados de batimetria e um polígono da Tanzânia e vejo o mesmo efeito. Veja minha atualização acima. Para contornos, é tão fácil como cont <- contour(r3)em seguida, lines(cont)não é de R grande?
Simbamangu
Ótimo, sim R pode ser fantástico! Obrigado @Simbamangu, isso funciona! No entanto, gostaria de traçar apenas meus polígonos espaciais com meu contorno de profundidade. Talvez crie um SpatialPolygonsDataframe? Usando intersect, recebo as informações de todo o quadrado que abrange meus polígonos espaciais. Eu preciso de alguma forma combinar meus polígonos espaciais com meu contorno profundo, mas não sei como fazê-lo. É por isso que pensei que a função 'extrair' era um bom começo, mas termino com uma lista de profundidade que não sei ao certo como posso me reconectar com meus polígonos espaciais ... Alguma idéia? Obrigado novamente!
Godina
Com o que você quer terminar? Contornos (SpatialLines) que se sobrepõem à sua subárea0? Você pode querer iniciar uma nova pergunta, pois agora isso está indo além do original.
Simbamangu #
Sim, eu gostaria de obter apenas minha subárea0 com o contorno de profundidade para esta área de definição: gis.stackexchange.com/questions/25112/…
GodinA
1

Aqui está outra solução, usando uma função básica de subconjunto em raster.

# create a raster with the dimensions that you are interested in.
r <- raster(extent(subarea0)+1) # +1 to increase the area
res(r) <- 0.1 # you can change the resolution here. 0.1 can be about 10x10 km if you are close to the equator.
crs(r) <- projection(subarea0) # transfer the coordinate system to the raster

# Now you can add data to the cells in your raster to mark the ones that fall within your polygon.
r[subarea0,] <- 1 # this an easy way to subset a raster using a SpatialPolygons object

# check your raster
plot(r)
ssanch
fonte