setwd(paste0(getwd(),"\\dados_curso"))
# carregar package:
library(raster)
library(sp)
library(rgdal) # para importar e exportar dados
library(rgeos) # para análise espacial
 
limite.ap<-as.matrix(read.table(file="limite.AP.SintraCascais.ETRS.txt",header=FALSE))
fich.mde<-"n38_w010_3arc_v2.tif"
mde<-raster(fich.mde)
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)
 
pov<-data.frame(
  designa=c("Sintra","Colares", "Alcabideche", "Almoçageme", "Azoia"),
  long=c(-9.3816589, -9.4426866999, -9.4098824999, -9.470629899, -9.476644970130),
  lat=c(38.8028687, 38.8066307, 38.7331569, 38.7967822, 38.7697973006),
  codPostal=c(2710,2705,2646,2705,2705)
)
 
wgs84<-"+proj=longlat +ellps=WGS84 +datum=WGS84"
xy<-cbind(pov$long,pov$lat) # matrix de coordenadas long/lat
 
xy.sp <- SpatialPoints(xy, proj4string=CRS(wgs84))
 
pov.wgs<- SpatialPointsDataFrame(xy.sp,pov[,c("designa","codPostal")])
 
coordinates(pov)<- ~ long + lat # converte em objecto sp
proj4string(pov) <- CRS(wgs84) # associa CRS
 
pov<-cbind(coordinates(pov),pov@data)
 
etrs<-"+proj=tmerc +lat_0=39.6682583 +lon_0=-8.1331083 +k=1 +x_0=0 +y_0=0 +ellps=GRS80 +units=m"
pov.etrs<-spTransform(pov.wgs,CRSobj=CRS(etrs))
 
coordinates(pov.etrs) # coordenadas dos pontos 
bbox(pov.etrs) #extensão (matrix 2*2)
 
pov.etrs@data # devolve data.frame
 
pov.etrs[2,] # devolve SpatialPointsDataFrame com uma única povoação (a 2a: Colares)
pov.etrs[pov.etrs$designa=="Colares",] # idem, mas usando atributo designa
 
pov.etrs[pov.etrs$codPostal>2700,] # devolve SpatialPointsDataFrame com 4 povoações
 
pov.etrs[,"designa"] # devolve SpatialPointsDataFrame com um único atributo (designa)
 
pov.etrs@data[,"designa"] # factor designa
 
as.character(pov.etrs@data[,"designa"]) # devolve vector do tipo caractere
 
# para definir as margens do plot
par(mfrow=c(1,1), mar=c(2,2,2,0))
plot(mde.etrs,  legend=FALSE, col=terrain.colors(100),xaxt="n", yaxt="n", box=FALSE, axes=FALSE, main="coordenadas ETRS PT TM 06")
# adicionar eixos com coordenadas
axis(1,at=round(bbox(pov.etrs)[1,]))
axis(2,at=round(bbox(pov.etrs)[2,]))
# adicionar polígono
polygon(limite.ap)
# escrever nomes de localidades
text(coordinates(pov.etrs), as.character(pov.etrs@data$designa),cex=.6)
 
part1 <- Polygon(limite.ap,hole=FALSE) 
 
mpart1 <- Polygons(list(part1),"SC") # criar objecto multiparte
 
ap.sp <- SpatialPolygons(list(mpart1),proj4string=CRS(etrs))
 
df <- data.frame(nome=c("Sintra Cascais"),row.names=c("SC"))
 
ap.etrs <- SpatialPolygonsDataFrame(ap.sp,df)
 
writeRaster(mde.etrs,file="mde_sintra_cascais.tif",overwrite=TRUE)
 
writeOGR(ap.etrs,dsn=getwd(),layer="Limite",driver="ESRI Shapefile",overwrite_layer=TRUE)
writeOGR(pov.etrs,dsn=getwd(),layer="Localidades",driver="ESRI Shapefile",overwrite_layer=TRUE)
 
writeOGR(pov.wgs,dsn="pov_sintra.kml",layer="localidades",driver="KML",overwrite_layer=TRUE)
 
library(maptools)
writeSpatialShape(pov.wgs,fn="Localidades")
 
icnf<-readOGR(dsn=getwd(),layer="AP_JUL_2014",encoding="ISO8859-1")
 
icnf@proj4string
 
length(icnf@polygons)
 
head(icnf@data)
 
indice<-grep(pattern="montejunto",icnf@data$NOME,ignore.case = TRUE) # o índice é 12
icnf@data[indice,] # verificar a escolha
 
plot(icnf[12,]) # icnf[12,] ainda é SpatialPolygonsDataFrame
 
icnf[12, ]@bbox # matrix com linhas "x" e "y" e colunas "min" e "max"
 
icnf@polygons[[12]]@ID  # é "11" porque os IDs começam em 0
 
icnf@polygons[[12]]@area 
 
icnf@polygons[[12]]@labpt 
 
icnf@data[12,"NOME"] # ou, equivalentemente,
icnf[12,]@data$NOME 
 
length(icnf@polygons[[12]]@Polygons) #tem 6 partes
 
for (i in 1:6) print(icnf@polygons[[12]]@Polygons[[i]]@hole)
 
for (i in 1:6) print(icnf@polygons[[12]]@Polygons[[i]]@area)
 
pol3 <- icnf@polygons[[12]]@Polygons[[3]]@coords # matriz com 2 colunas
 
# fazer figura  e colocar nome no canto xmin, ymax
plot(icnf[12,])
text(x=icnf[12,]@bbox["x","min"], 
      y=icnf[12,]@bbox["y","max"], 
      as.character(icnf@data[12,"NOME"]),pos=4,cex=.9) 
# representar a amarelo cada parte "negativa" e indicar área
for (i in 1:length(icnf@polygons[[12]]@Polygons)) 
{
     # aux é a i-ésima parte
     aux<-icnf@polygons[[12]]@Polygons[[i]];
     if (aux@hole) polygon(aux@coords,col="yellow")  
     text(x=aux@labpt[1], 
           y=aux@labpt[2],
           paste(round(aux@area/10000,1),"ha",sep=""),cex=.7,pos=4)
}
 
tipo.partes<-function(spdf)
{
  pos<-c(); neg<-c()
  for (i in 1:length(spdf@polygons)) # multipartes
  {
    p<-0; n<-0
    for (j in 1:length(spdf@polygons[[i]]@Polygons)) # partes
    {
      if (spdf@polygons[[i]]@Polygons[[j]]@hole) n<-n+1 else p<-p+1 # incrementa n se negativa; senão incrementa p
    }
    pos<-c(pos,p); neg<-c(neg,n)
  }
  return(list(positivos=pos,negativos=neg))
}
tipo.partes(icnf) # devolve resultado para dados do ICNF
 
library(ellipse)
 
multipartes<-list()
 
for (i in 1:length(icnf@polygons)) # i-ésima multiparte 
{
  partes<-list(); np<-0
  for (j in 1:length(icnf@polygons[[i]]@Polygons)) # j-ma parte 
  {
    if (icnf@polygons[[i]]@Polygons[[j]]@hole==FALSE)
    {
      np<-np+1
      partes[[np]] <- Polygon(round(ellipse(cov(icnf@polygons[[i]]@Polygons[[j]]@coords),centre = icnf@polygons[[i]]@Polygons[[j]]@labpt),2),hole=FALSE)
    }
  }
  multipartes[[i]]<-Polygons(partes, ID=icnf@polygons[[i]]@ID )
} 
 
elipses.etrs<-SpatialPolygons(multipartes,proj4string=icnf@proj4string)
elipses <- SpatialPolygonsDataFrame(elipses.etrs,data=icnf@data)
 
writeOGR(elipses,dsn=getwd(),layer="Elipses",driver="ESRI Shapefile",overwrite_layer=TRUE)
writeSpatialShape(elipses,fn="Elipses") # se houver dificuldades com writeOGR
 
elipses.wgs<-spTransform(elipses,CRS(wgs84))
writeOGR(elipses.wgs,dsn="elipses.kml",layer="areas_protegidas",driver="KML",overwrite_layer=TRUE)
 
concelhos<-readOGR(dsn=getwd(),layer="conc_1998",encoding="ISO8859-1") # cuidado com o CRS
concelhos@proj4string # verificar que não tem transformação de datum!!
 
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
concelhos.etrs<-spTransform(concelhos,CRS(etrs))
 
plot(concelhos.etrs)
plot(elipses,col="green",add=TRUE)
plot(icnf,col="red",add=TRUE)
 
parcelas<-readOGR(dsn=getwd(),layer="parcRibatejo",encoding="ISO8859-1")
parcelas@proj4string<-CRS(igeoe) # corrige CRS que estava incompleto
parcelas # para escrever sumário dos dados
head(parcelas) # para ver primeiras linhas da tabela de atributos
levels(parcelas$COD_AJUDA) # há 3 níveis "FOR", "OUT" e "SUP"
 
fichs <- list.files(pattern="banda[1-7]")
b<-brick(stack(as.list(fichs)))
 
par(mar=rep(0,4))
b.proj<-projectRaster(b,crs=parcelas@proj4string) # lento
plotRGB(b.proj,r=4,g=3,b=2,stretch="lin") 
plot(parcelas,col=c("green", "orange", "yellow")[parcelas$COD_AJUDA],add=TRUE)
legend(x=parcelas@bbox["x","min"],y=parcelas@bbox["y","max"],legend=levels(parcelas$COD_AJUDA),fill=c("green", "orange", "yellow")) 
 
ajudas<-data.frame(codigo=levels(parcelas$COD_AJUDA),designacao=c("Foragem", "Outros", "Superficies"))
 
tabela.merge.ajudas<-merge(x=parcelas@data,y=ajudas,by.x="COD_AJUDA",by.y="codigo",all.x=TRUE)
head(tabela.merge.ajudas)
 
parcelas.merge<-merge(x=parcelas,y=ajudas,by.x="COD_AJUDA",by.y="codigo",all.x=TRUE)
head(parcelas.merge)
 
par(mar=rep(0,4))
plotRGB(b.proj,r=4,g=3,b=2,stretch="lin") 
plot(parcelas.merge,col=c("green", "orange", "yellow")[parcelas.merge$COD_AJUDA],add=TRUE)
legend(x=parcelas.merge@bbox["x","min"],y=parcelas.merge@bbox["y","max"],legend=levels(parcelas.merge$designacao),fill=c("green", "orange", "yellow"),bg="white") 
 
expressao<-"Vinha de Uva de Mesa"
cultura<-parcelas[grep(pattern=expressao,parcelas$NOME_CULTU),]
 
plotRGB(b.proj,r=4,g=3,b=2,stretch="lin", ext=cultura@bbox) 
plot(cultura,add=TRUE,col="yellow")
 
out<-extract(b.proj,cultura[20,],weights=TRUE,normalizeWeights=FALSE,cellnumbers=TRUE)
outmat<-out[[1]] # 1a componente da lista de output
 
par(mar=rep(1,4))
plotRGB(b.proj,r=4,g=3,b=2,stretch="lin",ext=cultura[20,]@bbox)
plot(cultura[20,],col=rgb(r=1,b=0,g=1,alpha=.7),add=TRUE)
xy<-coordinates(b.proj)[outmat[,"cell"],]
text(xy[,1],xy[,2],outmat[,"weight"])
 
out<-extract(b.proj,cultura,weights=TRUE,normalizeWeights=FALSE) # lento
length(out) 
 
outmat<-out[[1]]
for (i in 2:length(out)) outmat<-rbind(outmat,out[[i]])
dim(outmat) # devolve o tamanho da matriz resultante
 
outcel <- outmat[outmat[,"weight"]>0.95,] 
head(outcel)
 
d<-function(xy0,xy){(xy0[1]-xy[1])^2+(xy0[2]-xy[2])^2}
 
vertice<-function(linha,pos)
{
  if (pos=="inicio") return(linha[1,])
  if (pos=="fim") return(linha[nrow(linha),])
}
 
empilha<-function(listapols) # empilha lista de "poli-linhas"
{
  # extremidade da primeira linha
  xy0<-vertice(listapols[[1]],pos="fim")
  mindist<- 10^10 # número muito grande
  for (i in 2:length(listapols)) # procura linha mais próxima 
  {
    xy1<-vertice(listapols[[i]],pos="inicio")
    xy2<-vertice(listapols[[i]],pos="fim")
    # inverte ordem da linha se necessário:
    if (d(xy0,xy2)<d(xy0,xy1))
    {
      listapols[[i]]<-listapols[[i]][nrow(listapols[[i]]):1,]
      xy1<-xy2
    }
    if (d(xy0,xy1)<mindist)  {index<-i; mindist<-d(xy0,xy1)}
  }
  listapols[[1]]<-rbind(listapols[[1]],listapols[[index]])
  listapols[[index]]<-NULL # elimina linha já usada
  return(listapols)
}
 
library(rgdal)
library(sp)
library(rgeos)
urlcaop<-"http://www.dgterritorio.pt/ficheiros/cadastro/caop/caop_download/conttrococaop2015_3"
download.file(url=urlcaop,destfile="caop2015.zip",mode="wb")
unzip(zipfile="caop2015.zip")
caop2015<-readOGR(dsn=getwd(),layer="Cont_Troco_CAOP2015",encoding="ISO8859-1")
freguesia<-"Alcântara"
# extrair features para essa freguesia e criar novo objecto SpatialLinesDataFrame
sldf<-caop2015[caop2015@data$Ent_Dir==freguesia |  caop2015@data$Ent_Esq==freguesia,]
 
# lista de poligonos (matrizes xy) que formam sldf:
listapols<-list()
# porque cada feature corresponde a uma única linha poligonal basta considerar a primeira linha de cada feature de sldf:
for (i in 1:length(sldf)) listapols[[i]]<-sldf@lines[[i]]@Lines[[1]]@coords
# empilhar linhas para criar poligono único
for (i in 1:(length(listapols)-1)) listapols<-empilha(listapols)
 
plot(sldf)
lines(listapols[[1]],col="red")
 
part1 <- Polygon(listapols[[1]],hole=FALSE) 
mpart1 <- Polygons(list(part1),freguesia) # criar objecto multiparte
freguesia.sp <- SpatialPolygons(list(mpart1),proj4string=sldf@proj4string)
df <- data.frame(nome=freguesia,row.names=freguesia)
freguesia.spdf <- SpatialPolygonsDataFrame(freguesia.sp,df)
# fazer figura do SpatialPolygonDataFrame criado com a indicação da área da freguesia:
plot(freguesia.spdf,main=freguesia)
xy<-coordinates(freguesia.spdf)
text(xy[1],xy[2],paste0(round(gArea(freguesia.spdf)/10^6,3),"km2")) 
 
