Skip to content

Instantly share code, notes, and snippets.

@alfcrisci
Created February 21, 2020 09:27
Show Gist options
  • Save alfcrisci/f295d8b4c5aea59c55a8f1c069c4acc5 to your computer and use it in GitHub Desktop.
Save alfcrisci/f295d8b4c5aea59c55a8f1c069c4acc5 to your computer and use it in GitHub Desktop.
process_UKCP12km=function(file,res=12000,bandN=1) {
require(raster)
require(proj4)
tmn<-suppressWarnings(raster(file,band=bandN))
lon=unique(coordinates(tmn)[,1])
lat=sort(unique(coordinates(tmn)[,2]))
xyGrid<-matrix(c(rep(lon,length(lat)),rep(lat,each=length(lon))),ncol=2)
xyGrid.wgs84<-proj4::project(xyGrid,proj="+proj=ob_tran +o_proj=longlat +lon_0=378 +o_lon_p=0 +o_lat_p=39.25 +a=6371229 +b=6371229 +to_meter=0.0174532925199 +wktext",inverse=T) #to project in latlong
xyGrid.wgs84<-SpatialPointsDataFrame(coords=xyGrid.wgs84,data=data.frame(z=values(flip(tmn,direction='y'))),proj4string=CRS('+init=epsg:4326'))
xyGrid_3035<-spTransform(xyGrid.wgs84,CRS('+init=epsg:3035'))
r_3035<-raster(nrows=round(diff(range(coordinates(xyGrid_3035)[,1]))/res),
ncols=round(diff(range(coordinates(xyGrid_3035)[,2]))/res)+1,
xmn=min(coordinates(xyGrid_3035)[,1]),
xmx=max(coordinates(xyGrid_3035)[,1]),
ymn=min(coordinates(xyGrid_3035)[,2]),
ymx=max(coordinates(xyGrid_3035)[,2]),
crs=CRS('+init=epsg:3035'),
vals=NA)
cells<-extract(r_3035,xyGrid_3035,cellnumber=T)[,1]
values(r_3035)[cells]<-values(flip(tmn,direction='y')) ### <<<--- the flip here is mandatory
r_final<-focal(r_3035,w=matrix(rep(1,9),ncol=3),fun=mean,pad=F,na.rm=T,NAonly=T)
return(r_final)
}
@alfcrisci
Copy link
Author

co-authors with maurizio marchi IBBR CNR

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment