[R-br] Novo problema para mudar projeção de um raster

ASANTOS alexandresantosbr em yahoo.com.br
Quinta Julho 25 23:32:07 BRT 2013


Boa noite Pessoal,

       Consegui resolver um problema e consegui outro, transformei as 
coordenadas latlong de um geotiff do Topodata em utm usando 
spTransform(). Como não consegui transformar diretamente um 
SpatialPixelsDataFrame do topodata apenas colocando o novo CRS e em 
vários posts não vi nenhuma solução para isso, tentei fazer da seguinte 
maneira, transformei apenas o @coords do topodata de geo para utm e 
depois tentei inserir com replace() os novos valores em utm dentro do 
SpatialPixelsDataFrame e ocorreu o seguinte erro:

Erro em replace(dem.sp em coords[1:10026, 1:2], dem.sp em coords, 
Utm em coords[1:10026,  :
   valores negativos não são permitidos na subscrição de matriz

     Não sei por onde continuar uma vez que minhas coordenadas geo são 
negativas e Mod[] não funcionou, segue CRM:

require(rgdal)
require(raster)

# Criar uma área de menor dimensão que a imagem inteira que abarque a 
região de interesse
xcc<-773759.1
ycc<-7841546
p.central<-cbind(xcc,ycc)

###Criando os vértices da área
coordV<-NULL
coordV <- 
rbind(coordV,cbind(p.central[,1]+c(-1500,1500,1500,-1500,-1500),p.central[,2]+c(1500,1500,-1500,-1500,1500)))
coordV
plot(coordV[,1],coordV[,2])
points(p.central[,1],p.central[,2], col="red")
#

# Cria um polígono com o contorno definido
bnds <- cbind(x=c(coordV[,1]), y=c(coordV[,2]))

# CRS UTM
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")))
proj4string(SP) = CRS("+proj=utm +zone=23 +south +datum=WGS84 +units=m 
+no_defs") ## Projeção

### CRS GEO
SPgeo<- spTransform(SP, CRS("+proj=longlat +datum=WGS84"))
#
## baixa e abre o objeto topodata
url=("http://www.dsr.inpe.br/topodata/data/geotiff/19S435SN.zip")
download.file(url, destfile = "19S435SN.zip")

### descompacta
system("unzip 19S435SN.zip")

# abrir o topodata declividade 19S435SN
dem <- raster('19S435SN.tif')
###

### atribui a projeção longlat  e datum WGS83 ao raster.
proj4string(dem) <- CRS("+proj=longlat +datum=WGS84")

#Cortar uma área menor de interesse
dem.crop <- crop(dem, extent(SPgeo), snap='out')
contorno.na <- setValues(dem.crop, NA)
contorno.r <- rasterize(SPgeo, contorno.na)
dem.mask <- mask(x=dem.crop, mask=SPgeo)
dem.sp <- as(dem.mask, "SpatialPixelsDataFrame")

## Mudar projeção de LatLong para Utm
LatLong <- data.frame(X =dem.sp em coords[,1], Y = dem.sp em coords[,2])
coordinates(LatLong) <- ~ X+Y
proj4string(LatLong) <- CRS("+proj=longlat +datum=WGS84")
Utm <- spTransform(LatLong, CRS("+init=epsg:29193"))


## Tentativa de substituir as coordenadas geo para utm dentro do objeto 
dem.sp
replace(dem.sp em coords[1:10026, 1:2], dem.sp em coords,Utm em coords[1:10026,1:2])

#
proj4string(dem.sp) <- CRS("+proj=utm +zone=23 +south +datum=WGS84 
+units=m +no_defs")
writeGDAL(dem.sp,"dem.tif")
#


Obrigado,

-- 
======================================================================
Alexandre dos Santos
Proteção Florestal
IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
Campus Cáceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
Cáceres - MT                      CEP: 78.200-000
Fone: (+55) 65 8132-8112 (TIM)   (+55) 65 9686-6970 (VIVO)
e-mails:alexandresantosbr em yahoo.com.br
         alexandre.santos em cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
======================================================================



Mais detalhes sobre a lista de discussão R-br