setwd(paste0(getwd(),"\\dados_curso"))
# carregar package:
library(raster)
library(rgdal)
# ler ficheiro:
pan <- raster("landsat8pan.tif")
cores<-gray(seq(0,1,length.out=100))
plot(pan,  xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,horizontal=TRUE)
 
pan # descreve o conjunto de dados 
pan@crs # devolve o sistema de coordenadas  em formato proj.4 que é +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
# exportar para ficheiro GeoTIFF
writeRaster(r, file="out.tif",format="GTiff",  overwrite=TRUE) 
 
# criar vector de 100 tons de cinzento entre preto (0) e branco (1)
cores<-gray(seq(0,1,length.out=100))
plot(pan,col=cores)
# tirar caixa e colocar a legenda na horizontal
plot(pan,  xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,horizontal=TRUE) 
 
v<-values(pan) # vector numérico
 
range(v)
 
pan@data@min
pan@data@max
 
out<-hist(v,main="valores vs intensidade")
fqs<-out$counts # frequências absolutas das classes
# representar a transformação linear de valores dos pixels em intensidade
lines(x=range(v),y=range(fqs),col="blue")
# representar os valores de intensidade
axis(4,at=seq(min(fqs),max(fqs),length.out=5),labels=seq(0,1,length.out=5))
mtext("intensidade", side=4,line=-1.5,col="blue")
 
plot(pan,  zlim=c(0,25000), xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,horizontal=TRUE) 
 
plot(pan,  zlim=c(0,quantile(v,.999)), xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,legend=FALSE) 
 
head(coordinates(pan))
res(pan)
extent(pan)
 
plot(pan, xlim=c(570000,575000), ylim=c(4315000,4325000), zlim=c(0,25000), xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,legend=FALSE) 
 
plot(pan, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)  
box <- drawExtent() # clicar em 2 cantos
# definir a extensão da imagem usando box 
plot(pan, ext=box, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)  
 
plot(pan, ext=box/2, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)  
 
plot(pan, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)  
zoom(pan, new=FALSE, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores) # clicar em 2 cantos
 
click(pan) # devolve valor de pixels; parar com <ESC>
 
pts <- locator() # devolve uma lista de coords; parar com <ESC>
polygon(cbind(pts$x,pts$y),border="red") # adiciona à imagem o polígono formado pelos pontos digitalizados
 
pan.crop <- crop(x=pan,y=pts) # a nova extensão é dada por pts
plot(pan.crop, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)  
 
igeoe <- "+proj=tmerc +lat_0=39.66666666666666  +towgs84=-283.1,-70.7,117.4,-1.16,0.06,-0.65,-4.1 +lon_0=1 +k=1 +x_0=200000 +y_0=300000 +ellps=intl  +pm=lisbon +units=m" # Coordenadas ``militares''
m.igeoe <- projectRaster(pan, crs=igeoe)
 
vazio <- raster(xmn=195000, xmx=205000, ymn=230000, ymx=245000,resolution=100,crs=igeoe)
m2 <- projectRaster(from=pan,to=vazio,method="ngb") 
m3 <- projectRaster(from=pan,to=vazio,method="bilinear")
 
par(mfrow=c(2,1))
plot(pan, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores,zlim=c(0,25000))
xy<-as.vector(pan@extent) 
# xpd para poder escrever fora da imagem
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy),pos=c(1,1,2,2),xpd=TRUE)
 
plot(m.igeoe, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores,zlim=c(0,25000))
xy<-as.vector(m@extent) # extensão xmin, xmax, ymin, ymax
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy),pos=c(1,1,2,2),xpd=TRUE)
 
igeoe.grid <- CRS("+proj=tmerc +lat_0=39.66666666666666  +nadgrids=ptLX_e89.gsb +lon_0=1 +k=1 +x_0=200000 +y_0=300000 +ellps=intl  +pm=lisbon +units=m")
m.grid <- projectRaster(pan, crs=igeoe.grid)
 
plot(m.grid, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores,zlim=c(0,25000))
xy<-as.vector(m.grid@extent)
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy),pos=c(1,1,2,2),xpd=TRUE)
 
# projectar raster em WGS84
wgs84<-"+proj=longlat +ellps=WGS84 +datum=WGS84"
w <- projectRaster(pan, crs=wgs84)
plot(w, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores,zlim=c(0,25000))
# Determinar a extensão
xy<-as.vector(w@extent) 
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy,4),pos=c(1,1,2,2),xpd=TRUE)
# Fazer zoom sobre parede da baragem
zoom(w, new=FALSE, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores) # clicar sobre 2 pontos
# Usar locator para obter coordenadas lat/long da parede da barragem
locator() # parar com <ESC>
 
limite.ap <- as.matrix(read.table(file="limite.AP.SintraCascais.ETRS.txt",header=FALSE))
 
mde <- raster("n38_w010_3arc_v2.tif")
 
etrs <- "+proj=tmerc +lat_0=39.6682583 +lon_0=-8.1331083 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m"
mde.etrs <- projectRaster(mde,crs=etrs)
mde.etrs <- crop(mde.etrs,limite.ap)
 
ndvi<-raster("ndviSintra.tif")
 
ndvi<-raster("ndviSintra.tif")
vazio<-raster(ext=mde.etrs@extent,crs=mde.etrs@crs,resolution=res(mde.etrs))
ndvi.etrs<-projectRaster(from=ndvi,to=vazio)
 
minx<-min(limite.ap[,1]); maxx<-max(limite.ap[,1])
miny<-min(limite.ap[,2]); maxy<-max(limite.ap[,2])
amostra2D<-cbind(sample(x=minx:maxx,size=2000),sample(x=miny:maxy,size=2000))
 
library(sp)
# vector TRUE/FALSE do tamanho de amostra:
dentro<-point.in.polygon(point.x=amostra2D[,1], point.y=amostra2D[,2],pol.x=limite.ap[,1],pol.y = limite.ap[,2])==1
# selecção dos pontos dentro de limite.ap:
amostra2D.ap<-amostra2D[dentro,]
 
y<-extract(ndvi.etrs,amostra2D.ap)
x<-extract(mde.etrs,amostra2D.ap)
 
par(mfrow=c(2,1),mar=rep(2,4))
plot(mde.etrs,  xaxt="n", yaxt="n", box=FALSE, axes=FALSE, col=cores,legend=FALSE)
polygon(limite.ap,border="red")
xy<-as.vector(mde.etrs@extent) 
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy,0),pos=c(1,1,2,2),xpd=TRUE,cex=0.8)
text(x=mean(xy[1:2]),y=xy[4],"mde.etrs",col="yellow",pos=1)
plot(ndvi.etrs, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE)
polygon(limite.ap,border="black")
xy<-as.vector(ndvi.etrs@extent) 
text(xy[c(1,2,1,1)],xy[c(3,3,3,4)],round(xy,0),pos=c(1,1,2,2),xpd=TRUE,cex=0.8)
text(x=mean(xy[1:2]),y=xy[4],"ndvi.etrs",col="black",pos=1)
points(amostra2D.ap,col="red")
 
plot(y~x,xlab="elevação (m)",ylab="ndvi",main="AP Sintra-Cascais",cex.lab=1.3)
xx<-log(x[x>0]) # tranformação logarítmica
yy<-log(y[x>0])
ajust<-lm(yy~xx) # ajustar regressão
curve(exp(ajust$coef[1])*x^ajust$coef[2],add=TRUE,col="red")
text(mean(range(x,na.rm=TRUE)),0,paste(" y = ", round(exp(ajust$coef[1]),3), " * x^",round(ajust$coef[2],3),sep=""))
 
fichs <- list.files(pattern="banda[1-7]")
s <- stack(as.list(fichs))
s # devolve sumário dos dados
b <- brick(s)
b # idem
nlayers(b) # devolve número de camadas em b
box <-extent(c(500000,  520000, 4310000, 4330000))
 
plotRGB(b,r=4,g=3,b=2,stretch="lin",ext=box) 
 
b@crs
 
e<-extent(as.vector(c( 513316, 514509, 4317568 ,4318417)))
landsat<-crop(b,e)
 
wgs84<-"+proj=longlat +ellps=WGS84 +datum=WGS84"
aux<-projectRaster(landsat, crs=wgs84)
# o slot @extent devolve um objecto de classe extent 
aux@extent # coordenadas lat/long do objecto landsat
zona.central<-aux@extent/2 # coordenadas lat/long da parte central da extensão dada por aux@extent 
 
install.packages("dismo")
library(dismo)
gm.wgs<-gmap(zona.central,type="satellite",lonlat=TRUE,rgb=TRUE,scale=2) 
 
gm.wgs
 
etrs<-"+proj=tmerc +lat_0=39.6682583 +lon_0=-8.1331083 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m"
gm.etrs<-projectRaster(gm.wgs,crs=etrs)
landsat.etrs<-projectRaster(landsat,crs=etrs)
 
par(mfrow=c(1,1))
# composição Google maps em cor verdadeira 
plotRGB(gm.etrs,r=3,g=2,b=1,stretch="lin") 
# composição Landsat RGB=432
# alpha define transparência
plotRGB(landsat.etrs,r=4,g=3,b=2,stretch="lin",ext=gm.etrs@extent,alpha=100,add=TRUE) 
 
ndvi <- (b$banda4 - b$banda3)/(b$banda4 + b$banda3)
 
ncell(ndvi[is.na(ndvi)])/ncell(ndvi) #devolve a proporção de pixels com NA; pode usar-se length em vez de ncell
ndvi[ndvi<0] <- 0 # substituição de valor
 
tabela <- cbind(c(-1,.2,.5),c(.2,.5,1),1:3)
ndvi.c <- reclassify(ndvi,rcl=tabela,right=TRUE)
 
min(b) # RasterLayer dos mínimos nas 6 bandas
range(b) # devolve RasterBrick com min e max 
 
cellStats(b,"mean") # devolve um vector com 6 componentes
 
median(b) #dá erro
b.mediana<-overlay(b,fun=median) # lento, mas funciona
 
urlzip<-"http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/N37W008.hgt.zip" # cobre SE de Portugal
 
download.file(url=urlzip,destfile="N37W008.hgt.zip",mode="wb") 
 
unzip(zipfile="N37W008.hgt.zip") 
 
srtm8<-raster("N37W008.hgt")
 
srtm <- merge(srtm8,srtm9)
srtm[srtm<=0]<-NA # para não representar elevações negativas
 
plot(srtm, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)
 
library(RColorBrewer) # para criar paletas
display.brewer.all() # ver possíveis paletas de cores
 
N<-10; pal<-"Spectral"; pie(rep(1,N), col=brewer.pal(n=N,pal))
 
plot(srtm, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, horizontal=TRUE, col=brewer.pal(n=11, "Spectral"))
 
points(y=37.095753,x=-8.683319,col="white")
text(y=37.095753,x=-8.683319,"Lagos",col="white",pos=3)
 
HTTP<-"http://dds.cr.usgs.gov/srtm/version2_1/SRTM3/Eurasia/"
for (tile in tiles)
{
  nomezip<-paste0(tile,".hgt.zip")
  download.file(url=paste0(HTTP,nomezip),destfile=nomezip,mode="wb") 
  unzip(zipfile=nomezip) 
  srtm<-raster(paste0(tile,".hgt"))
  if (tile==tiles[1]) mde<-srtm else mde<-merge(mde,srtm)
}
plot(mde)
 
inc<-terrain(srtm,opt="slope", units="radians", neighbors=8)
 
declives<-100*tan(values(inc)) # em %
range(declives,na.rm=TRUE) # devolve minimo e máximo
 
xyz<-cbind(coordinates(srtm),declives)
colnames(xyz)<-c("long","lat","declive")
# which() é necessário pois declives pode tomar valor NA
loc.max<-xyz[which(declives==max(declives,na.rm=TRUE)),]
 
plot(srtm, xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE, col=cores)
points(x=loc.max["long"],y=loc.max["lat"],col="red",cex=2)
 
# cria um objecto do tipo extent que é um rectangulo em lat/long em torno do ponto 
ext<-extent(c(loc.max[1]-0.003,loc.max[1]+0.003,loc.max[2]-0.003,loc.max[2]+0.003))
# gmap faz download da imagem Google
plot(gmap(x=ext,type="hybrid",lonlat=TRUE,rgb=FALSE,scale=2) )
 
inclin<-terrain(srtm,opt="slope", units="radians", neighbors=8)
orient<-terrain(srtm,opt="aspect", units="radians", neighbors=8)
ilumin<-hillShade(slope=inclin, aspect=orient, angle=40, direction=150)
 
extensao<-c(-9,-8.3,37.1,37.5)
plot(ilumin,ext=extensao,xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE,col=cores)
plot(srtm,ext=extensao,xaxt="n", yaxt="n", box=FALSE, axes=FALSE, legend=FALSE,col=cores,add=TRUE,alpha=.5)
 
srtm.crop<-crop(srtm,extent(c(-8.7,-8.4,37.1,37.4))) # Monchique
srtm.etrs<- projectRaster(srtm.crop,crs=etrs) # para obter coordenadas em metros
Sx <- focal(srtm.etrs,w=matrix(c(-1,-2,-1,0,0,0,1,2,1)/(8*res(srtm.etrs)[1]),ncol=3))
Sy <- focal(srtm.etrs,w=matrix(c(1,0,-1,2,0,-2,1,0,-1)/(8*res(srtm.etrs)[2]),ncol=3))
declive <- sqrt(Sx^2+Sy^2)
 
f <- matrix(1, nrow=3, ncol=3)
# índice de rugosidade do terreno 
TRI <- focal(srtm, w=f, fun=function(x, ...) sum(abs(x[-5]-x[5]))/8, pad=TRUE, padValue=NA)
# índice de posição topográfica 
TPI <- focal(srtm, w=f, fun=function(x, ...) x[5] - mean(x[-5]), pad=TRUE, padValue=NA)
# amplitude das diferenças de elevação
rough <- focal(srtm, w=f, fun=function(x, ...) max(x) - min(x), pad=TRUE, padValue=NA, na.rm=TRUE) 
 
# filtro não linear a aplicar à imagem:
moda<-function(m,na.rm=TRUE)
{
  freqs<-table(m)
  valores<-as.numeric(names(freqs))
  return(valores[which(max(freqs)==freqs)])
}
f<- matrix(1, nrow=3, ncol=3)
# focal é lento porque aplica a função moda pixel a pixel:
out <- focal(x=ndvi.etrs>.4, w=f, fun=moda, pad=TRUE, padValue=NA, na.rm=TRUE) 
 
pdf(file="ndvi-maior-do-40-filtro-moda.pdf",width=6,height=6)
par(mfrow=c(1,2)); 
plot(ndvi.etrs>.4,main="AP-SC NDVI>0.4",legend=FALSE); 
polygon(limite.ap,border="blue")
plot(out,main="NDVI>0.4 após filtro de moda",legend=FALSE); 
polygon(limite.ap,border="blue")
graphics.off()
 
