setwd(paste0(getwd(),"\\dados_curso")) # definir pasta de trabalho
library(raster)
library(sp)
library(rgdal) # para importar e exportar dados
library(rgeos) # para análise espacial
icnf<-readOGR(dsn=getwd(),layer="AP_JUL_2014",encoding="ISO8859-1")
 
valid<-gIsValid(icnf,byid=TRUE,reason=TRUE) 
which(valid!="Valid Geometry") #  13 55 58 61 inválidos
 
valid[13]  # "Ring Self-intersection[-95366.8066 -113617.4198]"
plot(icnf[13,]) # Arriba Fossil da Costa da Caparica
zoom(icnf[13,],ext=extent(c(-95370,-95360, -113620,-113610)),new=FALSE)
# para ver as coordenadas dos vértices:
xy<-icnf[13,]@polygons[[1]]@Polygons[[1]]@coords
text(xy,apply(round(xy),1,paste,collapse=","))
# para ver os indices dos vértices
zoom(icnf[13,],ext=extent(c(-95370,-95360, -113620,-113610)),new=FALSE)
text(xy,labels=1:nrow(xy),pos=c(1,2,3,4),cex=.8) # o problema é nos vértices 2, 3 e 4
 
icnf.novo<-icnf
part1 <- Polygon(xy[-c(2,3,4),],hole=FALSE) 
# icnf[13,] tem uma única multiparte que vamos subtituir por feature13:
feature13 <- Polygons(list(part1),row.names(icnf)[13]) # criar objecto multiparte
icnf.novo@polygons[[13]]<-feature13
zoom(icnf.novo[13,],ext=extent(c(-95370,-95360, -113620,-113610)),new=FALSE)
 
valid<-gIsValid(icnf.novo,byid=TRUE,reason=TRUE) 
which(valid!="Valid Geometry") # feature 13 já não é inválida
icnfv<-icnf.novo[valid=="Valid Geometry",] #só features válidas
 
gArea(icnfv) # devolve a área total
gArea(icnfv,byid=TRUE) # áreas da totalidade de ``features''
 
plot(icnfv,col="green",border=FALSE)
areas <- round(gArea(icnfv,byid=TRUE)/10000) # areas em ha
xy <- coordinates(gPointOnSurface(icnfv,byid=TRUE))
text(xy[areas > 20000,1],xy[areas>20000,2],paste(areas[areas>20000]), cex=0.7)
 
rios<-readOGR(dsn=getwd(),layer="rios",encoding="ISO8859-1")
concelhos<-readOGR(dsn=getwd(),layer="conc_2013",encoding="ISO8859-1") # cuidado com o CRS
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"
concelhos@proj4string<-CRS(igeoe) # corrige CRS que estava incompleto
valid<-gIsValid(concelhos,byid=TRUE,reason=TRUE)
which(valid!="Valid Geometry") # devolve índice 88
concelhos[88,] # Ílhavo
concelhosv<-concelhos[valid=="Valid Geometry",]
 
rios.etrs<-spTransform(rios,icnfv@proj4string)
concelhos.etrs<-spTransform(concelhosv,icnfv@proj4string)
 
rio<-rios.etrs[grepl(pattern="(R|r)io.*(M|m)ondego",rios.etrs@data$DESIGNACAO),]
 
its<-gIntersects(concelhos.etrs,rio,byid=TRUE)
# its é uma matriz lógica (TRUE/FALSE) 48*277
# 48 features para rio e 277 para concelhos.etrs
dim(its)
 
# definir vector de TRUE e FALSE para os concelhos 
aux<-as.logical(apply(its,2,max)) 
 
plot(concelhos.etrs, xlim=c(-95000, 151000), ylim=c(11000,136000))
plot(concelhos.etrs[aux,],add=TRUE, col="yellow")
xy<-coordinates(concelhos.etrs[aux,])
text(xy[,1],xy[,2],labels=concelhos.etrs@data[aux,"Municipio"],cex=.5,col="blue")
plot(rio,add=TRUE,col="blue")
 
parque<-icnfv[grepl(pattern='estrela',icnfv@data$NOME,ignore.case = TRUE),]
 
its<-gIntersects(concelhos.etrs,parque,byid=TRUE)
aux<-as.logical(apply(its,2,max)) # vector lógico
concelhos.etrs[aux,]
 
its<-gWithin(concelhos.etrs,parque,byid=TRUE)
aux<-as.logical(apply(its,2,max)) # vector lógico
concelhos.etrs[aux,]
 
plot(parque,col="green")
# concelhos que intersectam o Parque
its<-gIntersects(concelhos.etrs,parque,byid=TRUE)
aux<-as.logical(apply(its,2,max)) # vector logico
plot(concelhos.etrs[aux,],add=TRUE,border="red")
# escrever nomes dos concelhos seleccionados
xy<-coordinates(gPointOnSurface(concelhos.etrs[aux,],byid=TRUE))
text(xy[,1],xy[,2],concelhos.etrs@data[aux,"Municipio"],cex=0.7)
# concelhos totalmente contidos no parque:
its<-gWithin(concelhos.etrs,parque,byid=TRUE)
aux<-as.logical(apply(its,2,max)) # vector logico
plot(concelhos.etrs[aux,],add=TRUE,col="red")
# escrever nome do concelho seleccionado
xy<-coordinates(gPointOnSurface(concelhos.etrs[aux,],byid=TRUE))
text(xy[,1],xy[,2],concelhos.etrs@data[aux,"Municipio"],cex=0.7)
 
simplifica<-function(spdf)
{
  sp<-gCentroid(spdf,byid=TRUE) # devolve centroides das "features"
  return(SpatialPointsDataFrame(sp,spdf@data))
}
 
areas.ha <- round(as.vector(gArea(icnfv,byid=TRUE))/10000)
icnfg<-icnfv[areas.ha>1000,]
 
centroides.concelhos<-simplifica(concelhos.etrs)
d<-gDistance(centroides.concelhos,icnfg,byid=TRUE)
# atribuir nomes às linhas e colunas da matriz para facilitar a interpretação da mesma
colnames(d)<-as.character(centroides.concelhos$Municipio)
rownames(d)<-as.character(icnfg$NOME)
 
dc<-round(apply(d,2,min))
names(dc)<-as.character(centroides.concelhos$Municipio) 
 
dc[dc==max(dc)] # Mira, a 81 km
 
ind<-which(dc==max(dc)) # devolve o índice do concelho (coluna de d)
aux<-d[,ind] # distâncias desse concelho a todas as áreas protegidas
aux[aux==min(aux)] # devolve a área protegida à menor distância do concelho
 
plot(concelhos.etrs)
plot(concelhos.etrs[ind,],col="red",add=TRUE)
plot(icnfg,col="green",add=TRUE)
plot(icnfg[aux==min(aux),],col="blue",add=TRUE)
 
plot(gBuffer(icnfv[12,],byid=TRUE,width=100),border="red")
plot(icnfv[12,],add=TRUE)
 
library(RColorBrewer) # para criar paletas
display.brewer.all() # ver possíveis paletas de cores
 
par(mar=rep(0,4))
plot(gBuffer(icnfv[12,],width=1100),border=FALSE)
for (i in 10:1) plot(gBuffer(icnfv[12,],byid=TRUE,width=100*i),add=TRUE,col=brewer.pal(n=11, "Spectral")[i])
plot(icnfv[12,],add=TRUE,col="green")
 
distritos<-gUnaryUnion(concelhos.etrs, id = concelhos.etrs$Distrito)
 
names(distritos)
 
# determinar áreas dos concelhos com função area()
area(concelhos.etrs)
# substituir atributo area de concelhos por esses novos valores
concelhos.etrs@data$area<-area(concelhos.etrs)
# agregar as áreas dos concelhos do mesmo distrito:
aggregate(area~Distrito,data=concelhos.etrs,sum) 
 
df.distritos<-data.frame(aggregate(area~Distrito,data=concelhos.etrs,sum), row.names=names(distritos))
# Definir SpatialPolygonsDataFrame
distritos.etrs<-SpatialPolygonsDataFrame(distritos,data=df.distritos)
 
plot(distritos.etrs,col=brewer.pal(n=11, "Spectral"))
# area em ha de cada distrito
xy<-as.matrix(coordinates(distritos.etrs))
text(xy[,1],xy[,2],round(distritos.etrs$area/10000))
 
itsct<-gIntersection(parque,concelhos.etrs, byid=TRUE)
 
names(itsct) # "34 108" "34 132" "34 134" "34 135" "34 138" "34 170" "34 245"
 
spl<-strsplit(names(itsct),split=" ")
spl.matrix<-matrix(unlist(spl),ncol=2,byrow=TRUE)
 
df.icnf<-parque@data[spl.matrix[,1],c("NOME","CLASSIFICA")]
 
df.conc<-concelhos.etrs@data[spl.matrix[,2],c("Municipio","Distrito")]
 
df.itsct <- data.frame(cbind(df.icnf,df.conc),row.names=names(itsct))
itsct.etrs<-SpatialPolygonsDataFrame(itsct,data=df.itsct)
 
par(mar=rep(0,4))
plot(itsct.etrs,col=brewer.pal(n=length(itsct.etrs), "Spectral"))
# para as várias features:
for (i in 1:length(itsct.etrs))
{
  xy<-itsct.etrs@polygons[[i]]@labpt
  text(xy[1],xy[2],itsct.etrs@data[i,"Municipio"],pos=3,cex=0.6)
  text(xy[1],xy[2],itsct.etrs@data[i,"Distrito"],pos=1,cex=0.6)
}
 
itsct.etrs@data # para observar a tabela de atributos
int<-itsct.etrs[itsct.etrs$Municipio=="OLIVEIRA DO HOSPITAL"  & itsct.etrs$NOME=="Serra da Estrela",]
gArea(int) # menos de 20 ha
 
novoParque<-gDifference(parque,int,byid=TRUE)
plot(parque,xlim=c(27142, 33728), ylim=c(73244,68468))
plot(novoParque,col="red",add=TRUE)
 
plot(icnfv,col="green",border=FALSE)
plot(gConvexHull(icnfv,byid=TRUE),add=TRUE)
 
library(rgdal)
# ler ficheiro
g <- readGDAL("landsat8pan.tif")
 
summary(g) # descreve g
str(g) # devolve a estrutura de g
 
# g@data é uma das "slots" do objecto: é uma data.frame
# os nomes das colunas de g@data são então dados por
names(g@data) # band1
# histograma dos valores da imagem:
hist(g@data$band1)
cores <- gray(seq(0,1,length.out=100)) 
# image é melhor do que plot, que bloqueia para landsat8 pan (225*10^6 pixels) 
image(g,col=cores,zlim=quantile(g@data$band1,probs=c(0,.99)))
 
writeGDAL(g,fname="out.tif") # 
 
rasterize(x=icnf[6,],y=mde.etrs) # icnf[6,] é a AP Sintra-Cascais
 
ogrDrivers()
 
writeFormats()
 
