spatial R

exampleSpatialPolygons.R

library(rgeos)
library(maptools)
trentino_crs=CRS('+init=epsg:32632')
#create a list of polygomns
spatialPolygonList=lapply(seq_along(wkt_list), function(i) {
    a=readWKT(wkt_list[i,]$WKT,p4s=trentino_crs)
    a@polygons[[1]]@ID=as.character(i)
   return(a)
})
# join polygons in 1 SpatialPolygons
joined = SpatialPolygons(lapply(spatialPolygonList, function(x){x@polygons[[1]]}))
# create a spatial polygon from extent
myPoly<-gBuffer(bbox2SP(bbox=bbox(joined),proj4string=trentino_crs), width=buffer_width, byid=TRUE)

extractSlopeAndAspectValues = function(dem,mask){  
    execGRASS("g.region", parameters=list(raster=mask,res='2'))
    execGRASS("r.mask", parameters=list(raster=mask),flags=c("overwrite"))
    elevation = raster(readRAST(c(dem)))
    print("calcolo aspect")
    x = terrain(elevation, opt=c('aspect'), unit='degrees', neighbors=8)
    esposizione_value=getValues(x$aspect)
    esposizione = getAspectFromNorth(mean(esposizione_value,na.rm=T))
    esposizione_005_quantile=getAspectFromNorth(quantile(esposizione_value,0.05,na.rm=T))
    esposizione_095_quantile=getAspectFromNorth(quantile(esposizione_value,0.95,na.rm=T))
    x = terrain(elevation, opt=c('slope'), unit='tangent', neighbors=8)
    slope_value=getValues(x$slope)
    slope = mean(slope_value,na.rm=T)
    slope_005_quantile=quantile(slope_value,0.05,na.rm=T)
    slope_095_quantile=quantile(slope_value,0.95,na.rm=T)
    return(c(slope,slope_005_quantile,slope_095_quantile,esposizione,esposizione_005_quantile,esposizione_095_quantile));
}



getAstroSunData = function (lon,lat,year){
  start_date <- as.POSIXct(paste(year,"-01-01", tz="Europe/Rome"))
  sequence <- seq(from=start_date, length.out=365 , by="days")
  sunrise <- sunriset(lon.lat, sequence, direction="sunrise", POSIXct.out=TRUE)
  sunset <- sunriset(lon.lat, sequence, direction="sunset", POSIXct.out=TRUE)
  data.frame(date=as.Date(sunrise$time),
             sunrise=as.numeric(format(sunrise$time, "%H%M")),
             sunset=as.numeric(format(sunset$time, "%H%M")),
             day_length=as.numeric(sunset$time-sunrise$time))

            
    }