Extrair a varredura da varredura usando o shapefile do polígono em R

13

Eu sou novo no R e estou usando o pacote raster. Tenho um problema ao extrair polígonos de um arquivo raster existente. Se eu usar

extract(raster, poly_shape)

Na varredura, sempre cria uma lista com os dados. O que eu realmente quero é extrair outro arquivo raster que eu possa carregar com o ArcGIS novamente. Depois de ler um pouco mais, acho que a função de recorte é o que realmente preciso. Mas quando eu tento usar essa função

crop(raster, poly_shape)

Eu recebo este erro:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

Os arquivos raster e poly_shape são iguais para as duas funções. Você pode me dizer o que poderia estar errado aqui? É certo que a função de corte crie outra varredura e não uma lista?

Edição : A função extensão () não funciona para mim. Ainda obtenho o mesmo erro. Mas tenho certeza que os 2 conjuntos de dados se sobrepõem! Com o

extract(raster, poly_shape)

Eu recebo os dados certos. Apenas como uma lista e não como uma varredura como eu quero tê-la. Acabei de carregar os conjuntos de dados no ArcGIS antes e eles se encaixam muito bem, portanto não verifiquei a projeção. Agora eu tentei

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

e você pode ver que as projeções não se encaixam. A função extrair parece capaz de transformar automaticamente os arquivos da maneira correta. Eu sei disso porque fiz o seguinte:

  • Eu cortei a parte exata do polígono que extraí no R também no ArcGIS
  • Calculei a soma de todos os valores do polígono R extraído (lista)
  • Calculei a soma de todas as células raster que recortei no ArcGIS

Os 2 têm exatamente o mesmo resultado, então acho que a conclusão deve ser que a função extrair funcionou corretamente. Agora eu tenho duas opções, eu acho:

  1. Preciso de uma maneira de tirar um Raster da lista extraída novamente ou
  2. Os 2 conjuntos de dados (raster + poly_shape) precisam usar a mesma projeção e a função de corte deve funcionar

O que você sugeriria fazer aqui?

Lars
fonte
E se for um rgbi de 4 bandas? As bandas estão perdidas até agora ...
Doris
Bem-vindo ao GIS SE! Como um novo usuário, faça um breve tour . Em seguida, considere editar sua resposta para fornecer informações e referências adicionais. Consulte Como escrevo uma boa resposta? para mais informações.
Andy
Se você tiver uma nova pergunta, faça-o clicando no botão Fazer pergunta . Inclua um link para esta pergunta se ela ajudar a fornecer contexto. - Do comentário
Erik

Respostas:

19

A função extrair está se comportando exatamente como deveria. Você pode forçar a função de corte a usar a extensão do polígono e depois mascarar o objeto para retornar a varredura exata que representa a área do polígono. Se você continuar recebendo o erro, isso significa que seus dados não se sobrepõem.

Lembre-se de que o R não realiza projeção "on the fly"; portanto, verifique suas projeções. Você pode verificar se suas extensões se sobrepõem usando a função "extension ()".

Aqui está um exemplo de corte usando a extensão do polígono e mascarando a varredura resultante usando o polígono "rasterizado".

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)
Jeffrey Evans
fonte
4
extract () não executar no reprojeção mosca, mas de culturas () não. Essa conta pode por alguma confusão
mdsumner
@Jefferey crop () e mask () apenas cortam a varredura de acordo com as extensões retangulares do polígono que não a cortam dentro dos limites do polígono. Alguma idéia de quais comandos podem cortar a varredura dentro dos limites do polígono fornecido?
Csheth
1
@Chintan Sheth, para que a máscara subconjunto dentro do polígono, você precisa ter uma varredura representando os valores dentro do polígono. É por isso que você rasteriza o polígono do subconjunto e depois o mascara. A etapa de corte é reduzir a extensão da varredura para a mesma do polígono do subconjunto, que no passado, se ignorado, geraria um erro de incompatibilidade de extensão.
Jeffrey Evans
spTransformdo sppacote (que às vezes é carregado automaticamente com outros pacotes R espaciais) permite reprojetar para que os dois arquivos estejam na mesma projeção, por exemplo. good_poly=spTransform(spolygon, CRSobj=crs(raster_file))
User3386170
@ user3386170, hein? Não sei ao que você está chegando. Essa pergunta ocorreu no momento em que o pacote raster acabou de adicionar "projeção instantânea" em algumas funções. Essas funções geraram um erro anteriormente quando as projeções não correspondiam, mas esta postagem foi de 2014. Você também deve estar sempre ciente de sempre carregar o rgdal ao usar sp :: spTransform (), pois adiciona funcionalidade adicional importante aos bastidores.
9118 Jeffrey Evans
8

O que eu realmente procurei foi a mask()função.

mask(raster, poly_shape)

funciona sem erros e retorna o que procurei.

Lars
fonte
2
Reprojete seus dados para que estejam no mesmo espaço de projeção. Mesmo no ArcGIS, onde a projeção instantânea é automática, é uma prática muito ruim realizar análises em diferentes projeções. Com dados em projeções diferentes, seus dados não compartilharão extensões comuns, que é o erro que você está recebendo.
11134 Jeffrey Evans
Use projectExtent () para obter apenas a extensão para cortar a varredura.
Mdsumner #
Para ajustar o formato de perguntas e respostas do site, isso deve ser colocado no corpo principal da pergunta como uma edição / atualização (e adicione um comentário à resposta, em "resposta", para que eles saibam que há mais coisas a serem examinadas).
Matt Wilkie
@mattwilkie Desculpe por não ajustar o formato, mas meu texto foi muito longo para publicá-lo como um comentário aqui. @JeffreyEvans Eu tentei o seguinte: projection(raster) = projection(poly_shape)e o contrário projection(poly_shape) = projection(raster), mas em ambos os sentidos produzir o mesmo erro: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. Existe uma maneira de ver qual projeção é usada em tempo real usando a função extract () (porque essa obviamente funciona)?
Lars
1
O que eu realmente procurei foi a função mask (). mask(raster, poly_shape)funciona sem erros e retorna o que procurei.
Lars
-1

A extensão funciona muito bem ... Eu acho que o Xmin, Xmax, Ymin e Ymax da sua extensão são diferentes dos X e Y da sua varredura - ou seja, eles estão no lado oposto

ToNoY
fonte