knitr::opts_chunk$set(echo = TRUE) require(hydroCR) require(rcm2019) fill.na <- function(x, i=5) { if( is.na(x)[i] ) { return( round(mean(x, na.rm = TRUE), 4) ) } else { return( x[i] ) } } preciwei = function(x){ w = x/sum(x) return(w/sum(w)) } ################ setwd('/home/mha/Dropbox/melichar/data/GIS/') s = readOGR('PLO40hranice_S-JTSK_Krovak_East_North.shp') ws = krov2wgs(s) rws = rotatePoly(ws) dem = raster('srtm_krov500r.tif') projection(dem) = s@proj4string asp = terrain(dem, opt = 'aspect') slo = terrain(dem, opt = 'slope') xy = coordinates(dem) xy = krov2wgs(SpatialPoints(xy)) rxy = rotatePoints(xy) setdp('chars') preci = raster('srazky_8110.nc') kpreci = projectRaster(preci, crs = s@proj4string) #temp = raster('teplota_8110.tif') #ktemp = projectRaster(temp, crs = s@proj4string) preci = crop(preci, buffer(ws, 0.2)) pxy = coordinates(preci) rpxy = rotatePoints(SpatialPoints(pxy)) rpxy = SpatialPointsDataFrame(rpxy, data = data.frame(PR = values(preci))) finepreci = resample(kpreci, dem, method = 'ngb') #ktemp = crop(ktemp, buffer(s, 10000)) setwd(file.path(Sys.getenv('R_DATA_PATH'), 'RCM2019', 'annual')) #setwd('/media/mha/KINGSTON/RCM2019/') d = dir(pattern = 'tas_') f = d[1] setwd(file.path(Sys.getenv('R_DATA_PATH'), 'RCM2019', 'annual')) p = brick(f, varname = 'pr')
Původní RCM data (viz obr. níže vlevo) byla bilineárně interpolována do 1 km rozlišení (obr. vpravo)
par(mfrow = c(1, 2)) plot(p[[1]] ) plot(disaggregate(p[[1]], fact = 10, method = 'bilinear'))
Následně byla vytvořena na základě dem korekční vrstva opravující teplotu o 0.65 st. C / 100 m nadm. výšky (níže vlevo). Tato vrstva byla použita ke korekci RCM interpolované vrstvy.
aele = cellStats(mask(dem, s), mean) cdem = ((aele-dem)/100)*0.65 par(mfrow = c(1, 2)) plot(cdem) plot(disaggregate(p[[1]], fact = 10, method = 'bilinear'))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.