Re: [R-br] Fwd: como extrair dados de arquivos GEOtiff

Olá, boa tarde! Sugiro tabalhar com a biblioteca {rgdal} (e {sp}), pois ela já interpreta o cabeçalho geotiff, lendo seus dados na escala apropriada. Segue uma rotina que usei pra ler dados de altitude tomados do projeto TOPODATA (INPE). Acredito que possa ser facilmente adaptada ao teu caso. Preciso dar uma melhorada, mas acho que dá pra pegar a ideia. ### <begin> ### #install.packages(c("geoR", "maptools","sp", "rgdal"), dep=TRUE) #require(geoR); require(maptools); require(sp); require(rgdal) ### Download do exemplo (399 Kb) para a pasta definida em getwd() download.file('https://dl.dropbox.com/u/117618178/dourados/DDOS.tif', file.path(getwd(),'DDOS.tif'), mode='wb') alt = readGDAL("DDOS.tif"); gridded(alt) ### 'importa' o geotiff image(alt) ### visualiza na forma de imagem image(alt, col=terrain.colors(21)) ## outra opção de visualização slotNames(alt) ### nome dos slots (S4) proj4string(alt) ### para conferir a projeção (definida no geotiff) head(alt@data) ### observe os valores dos dados (neste caso, altitudes) ### se quiser todos os valores, alojando em um objeto DF alt.sp <- SpatialPoints(alt); slotNames(alt.sp) alt.df <- data.frame(x=alt.sp@coords[,1],y=alt.sp@coords[,2],alt=alt@data) head(alt.df); length(alt.df$x) ### extrair as altitudes em alguns pontos aleatórios (ou não! ;D) ### exemplo com 14 pontos pts <- structure(list(x = c(712930.163349833, 715345.025457117, 718363.603091221, 720636.414486312, 724152.169613092, 723903.580866754, 721701.79482776, 719215.90736438, 716126.304374179, 713391.828164461, 713924.518335185, 717582.324174159, 720352.313061925, 723086.789271644), y = c(7548164.89958582, 7548306.95029801, 7548271.43761996, 7548271.43761996, 7548129.38690777, 7546105.16425902, 7546105.16425902, 7546069.65158097, 7546211.70229316, 7546034.13890292, 7543441.71340539, 7543583.76411759, 7543654.78947368, 7543832.35286393)), .Names = c("x", "y")) ### Coloca tudo na classe Spatial pra poder manipular com funções do {sp} pts.sp <- SpatialPoints(pts) ### transforma os pontos para SpatialPoints pts.sp@proj4string;alt@proj4string ### observe as projeções definidas pts.sp@proj4string <- alt@proj4string ### atribui projeção dos pontos usando do grid over(pts.sp, alt) ### extrai valor em cada par de coordenadas ### extrair as altitudes usando uma grade regular ### neste caso vou ler a cada 300m alt@bbox ### limites da imagem/grade de altitude grd <- expand.grid(seq(711985,725005,300), seq(7541995,7549015,300)) ### modo 1 grd <- expand.grid(seq(alt@bbox[1,1],alt@bbox[1,2],300), seq(alt@bbox [2,1],alt@bbox[2,2],300)) ### modo 2 grd.sp <- SpatialPoints(grd) grd.sp@proj4string;alt@proj4string ### observe as projeções definidas grd.sp@proj4string <- alt@proj4string ### atribui projeção dos pontos usando do grid over(grd.sp, alt) ### xtrai valor em cada par de coordenadas ### <end> ### ================================================ Éder Comunello Ph.D. Student in Agricultural Systems Engineering (USP/ESALQ) Piracicaba, SP, Brazil [22 42.7'S, 47 37.8'W] Researcher at Embrapa Western Region Agriculture Dourados, MS, Brazil [22 16.5'S, 54 49.0'W] ================================================ UTC-03:00

Olá Éder, Sua rotina funcionou perfeitamente aqui com meus dados!!! Wow, isso foi uma mão na roda, muitíssimo obrigado! Att Carlos Em 14 de novembro de 2012 18:33, Eder Comunello <ecomunel@gmail.com>escreveu:
Olá, boa tarde!
Sugiro tabalhar com a biblioteca {rgdal} (e {sp}), pois ela já interpreta o cabeçalho geotiff, lendo seus dados na escala apropriada.
Segue uma rotina que usei pra ler dados de altitude tomados do projeto TOPODATA (INPE). Acredito que possa ser facilmente adaptada ao teu caso.
Preciso dar uma melhorada, mas acho que dá pra pegar a ideia.
### <begin> ### #install.packages(c("geoR", "maptools","sp", "rgdal"), dep=TRUE) #require(geoR); require(maptools); require(sp); require(rgdal)
### Download do exemplo (399 Kb) para a pasta definida em getwd() download.file('https://dl.dropbox.com/u/117618178/dourados/DDOS.tif', file.path(getwd(),'DDOS.tif'), mode='wb')
alt = readGDAL("DDOS.tif"); gridded(alt) ### 'importa' o geotiff image(alt) ### visualiza na forma de imagem image(alt, col=terrain.colors(21)) ## outra opção de visualização slotNames(alt) ### nome dos slots (S4) proj4string(alt) ### para conferir a projeção (definida no geotiff) head(alt@data) ### observe os valores dos dados (neste caso, altitudes)
### se quiser todos os valores, alojando em um objeto DF alt.sp <- SpatialPoints(alt); slotNames(alt.sp) alt.df <- data.frame(x=alt.sp@coords[,1],y=alt.sp@coords[,2],alt=alt@data) head(alt.df); length(alt.df$x)
### extrair as altitudes em alguns pontos aleatórios (ou não! ;D) ### exemplo com 14 pontos pts <- structure(list(x = c(712930.163349833, 715345.025457117, 718363.603091221, 720636.414486312, 724152.169613092, 723903.580866754, 721701.79482776, 719215.90736438, 716126.304374179, 713391.828164461, 713924.518335185, 717582.324174159, 720352.313061925, 723086.789271644), y = c(7548164.89958582, 7548306.95029801, 7548271.43761996, 7548271.43761996, 7548129.38690777, 7546105.16425902, 7546105.16425902, 7546069.65158097, 7546211.70229316, 7546034.13890292, 7543441.71340539, 7543583.76411759, 7543654.78947368, 7543832.35286393)), .Names = c("x", "y"))
### Coloca tudo na classe Spatial pra poder manipular com funções do {sp} pts.sp <- SpatialPoints(pts) ### transforma os pontos para SpatialPoints pts.sp@proj4string;alt@proj4string ### observe as projeções definidas pts.sp@proj4string <- alt@proj4string ### atribui projeção dos pontos usando do grid over(pts.sp, alt) ### extrai valor em cada par de coordenadas
### extrair as altitudes usando uma grade regular ### neste caso vou ler a cada 300m alt@bbox ### limites da imagem/grade de altitude grd <- expand.grid(seq(711985,725005,300), seq(7541995,7549015,300)) ### modo 1 grd <- expand.grid(seq(alt@bbox[1,1],alt@bbox[1,2],300), seq(alt@bbox [2,1],alt@bbox[2,2],300)) ### modo 2 grd.sp <- SpatialPoints(grd) grd.sp@proj4string;alt@proj4string ### observe as projeções definidas grd.sp@proj4string <- alt@proj4string ### atribui projeção dos pontos usando do grid over(grd.sp, alt) ### xtrai valor em cada par de coordenadas ### <end> ###
================================================ Éder Comunello
Ph.D. Student in Agricultural Systems Engineering (USP/ESALQ) Piracicaba, SP, Brazil [22 42.7'S, 47 37.8'W]
Researcher at Embrapa Western Region Agriculture Dourados, MS, Brazil [22 16.5'S, 54 49.0'W] ================================================ UTC-03:00
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- _____________________________________ *Carlos D'Apolito - PhD candidate - ID:1239480* Earth Sciences School of Geography & Environmental Sciences Aston Webb Building University of Birmingham Birmingham B15 2TT UK phone: 44 - 07 919 673 173 CXD280@bham.ac.uk CV: http://lattes.cnpq.br/4557754746424026
participantes (2)
-
D'Apolito
-
Eder Comunello