Obrigado Mauro, Resolvi o erro na rotina e mesmo assim o spTransform() não esta convertendo de utm na zona 23K para latlong no hemisfério sul, segue CRM completo ate o recorte da área de interesse no geoTiff do topodata: require(sp) 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))) 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 em lat long para o polígomno ficar na mesma projeção que o GeoTif SPlatlong<- spTransform(SP, CRS("+proj=longlat +south +datum=WGS84")) # # #Obtenção do objeto declividade do TOPODATA url=("http://www.dsr.inpe.br/topodata/data/geotiff/19S42_SN.zip") download.file(url, destfile = "19S42_SN.zip") ### descompacta system("unzip 19S42_SN.zip") demlatlog<-raster('19S42_SN.tif')### Abre o slope TOPODATA # #Cortar a região de interesse demlatlog.crop <- crop(demlatlog, extent(SPlatlong), snap='out') # Obrigado, Alexandre Em 22/07/2013 22:09, Mauro Sznelwar escreveu:
Tentei rodar o seu CRM e deu problema aqui
bnds <- cbind(x=c(coordV2[21:25,1]), y=c(coordV2[21:25,2])) Error in coordV2[21:25, 1] : subscript out of bounds
Segue CRM:
require(sp) require(rgdal) # 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 <- rbind(coordV,cbind(p.central[,1]+c(-1500,1500,1500,-1500,-1500),p.central[,2]+c(1500,1500,-1500,-1500,1500))) coordV coordV2=coordV[!is.na(coordV[,1]),] #plot(coordV2[21:25,1],coordV2[21:25,2]) #points(p.central[,1],p.central[,2], col="red") #
# Cria um polígono com o contorno definido bnds <- cbind(x=c(coordV2[21:25,1]), y=c(coordV2[21:25,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 em lat long para o polígomno ficar na mesma projeção que o GeoTif SPlatlong<- spTransform(SP, CRS("+proj=longlat +south +datum=WGS84")) #
Obrigado,
Alexandre
-- ====================================================================== 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@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================