Generating the 5 km RHEAS setup - SERVIR/RHEAS GitHub Wiki
Let’s set the bounding box for our domain
minlat = 5.5
maxlat =28.75
minlon = 92.0
maxlon = 109.75Now let’s download the land cover dataset from MODIS and the GTOPO30 elevation dataset from the USGS Earth Explorer, from which we will get a set of HDF and TIF files. We need to reproject, mosaic and subset them before we continue
for f in *hdf; do gdalwarp -t_srs epsg:4326 HDF4_EOS:EOS_GRID:"$f":MCD12Q1:LC_Type2 `echo $f | sed 's/hdf/tif/'`;done
gdal_merge.py -o lc.tif MCD12Q1*.tif
gdal_translate -projwin 91.0 29.75 110.75 4.5 lc.tif lmr_lc.tif
gdal_merge.py -o gt30.tif gt30*.tif
gdal_translate -projwin 91.0 29.75 110.75 4.5 gt30.tif lmr_dem.tifThen we enter the scripts directory of RHEAS and generate the 5-km files
from resample_vic import *
resample_soil(5.5, 28.75, 92.0, 109.75, 0.05, "../../global_soil_0.25deg.txt")
resample_vegetation("../../lmr_lc.tif", "../../global_soil_0.05deg.txt", "../../vic_veglib.txt", "../../global_lai_0.05deg.txt", 0.05)
resample_snowbands("../../lmr_dem.tif", "../../global_soil_0.05deg.txt", "../../global_snowbands_0.05deg.txt", 0.05, 5)The next steps involve updating the database so it knows about the new files. First, we insert the file information
echo "insert into vic.input values (0.05, 'vic/global_snowbands_0.05deg.txt', 'vic/global_lai_0.05deg.txt', 'vic/vic_veglib.txt', 'vic/global_soil_0.05deg.txt', 2, 'vic/dssat.inp.base')" | ./bin/psql -d rheasand then we insert each grid cell’s information into the soils database table
echo "begin;" > vicsoils5km.sql
awk '{printf("insert into vic.soils (id,line,depths,elev,geom,resolution) values (%d, '\''%s'\'', ARRAY[%f,%f,%f],%.2f,st_geomfromtext('\''POINT(%f %f)'\'',4326),0.05)\n;",$2,$0,$23,$24,$25,$22,$4,$3)}' ../data/vic/global_soil_0.05deg.txt >> vicsoils5km.sql
echo "end;" >> vicsoils5km.sql
./bin/psql -d rheas -f scripts/vicsoils5km.sqlThe last step involves creating the cached 5-km raster tables for the meteorological datasets
echo "create table precip.chirps_20 as (with f as (select fdate,st_rescale(st_tile(rast,10,10),0.05) as rast from precip.chirps) select fdate,rast,dense_rank() over (order by fdate,st_upperleftx(rast),st_upperlefty(rast)) as rid from f);" | ./bin/psql -d rheas
echo "create table tmax.ncep_20 as (with f as (select fdate,st_rescale(st_tile(rast,10,10),0.05,'bilinear') as rast from tmax.ncep) select fdate,rast,dense_rank() over (order by fdate,st_upperleftx(rast),st_upperlefty(rast)) as rid from f);" | ./bin/psql -d rheas
echo "create table tmin.ncep_20 as (with f as (select fdate,st_rescale(st_tile(rast,10,10),0.05,'bilinear') as rast from tmin.ncep) select fdate,rast,dense_rank() over (order by fdate,st_upperleftx(rast),st_upperlefty(rast)) as rid from f);" | ./bin/psql -d rheas
echo "create table wind.ncep_20 as (with f as (select fdate,st_rescale(st_tile(rast,10,10),0.05,'bilinear') as rast from wind.ncep) select fdate,rast,dense_rank() over (order by fdate,st_upperleftx(rast),st_upperlefty(rast)) as rid from f);" | ./bin/psql -d rheasassuming that we want to only resample the CHIRPS and NCEP datasets.