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

[Eder Comunello]

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 em 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 em coords[,1],y=alt.sp em coords[,2],alt=alt em 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 em proj4string;alt em proj4string ### observe as projeções definidas
pts.sp em proj4string <- alt em 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 em 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 em bbox[1,1],alt em bbox[1,2],300), seq(alt em bbox
[2,1],alt em bbox[2,2],300)) ### modo 2
grd.sp <- SpatialPoints(grd)
grd.sp em proj4string;alt em proj4string ### observe as projeções definidas
grd.sp em proj4string <- alt em 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
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121114/c0643150/attachment.html>