forked from edwinkost/extreme_value_analysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
making_landmask_for_extreme_value_analysis.sh
40 lines (31 loc) · 2.03 KB
/
making_landmask_for_extreme_value_analysis.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
set -x
# making the landmask map for extreme value analysis
cd /scratch-shared/edwinsut/clone_for_extreme_value_analysis/
rm -r landmask_extreme_value_analysis
mkdir landmask_extreme_value_analysis
cd landmask_extreme_value_analysis
# identify cells with forcing
gdal_translate -of PCRaster ../forcing_extent/pr_gfdl_1_jan_2006.nc pr_gfdl_1_jan_2006.map
gdal_translate -of PCRaster ../forcing_extent/pr_watch_1_jan_1958.nc pr_watch_1_jan_1958.map
pcrcalc landmask_forcing.map = "if(if(defined(pr_gfdl_1_jan_2006.map), defined(pr_watch_1_jan_1958.map)), boolean(1.0))"
mapattr -p *.map
aguila landmask_forcing.map
# resample landmask of forcing to 5 arcmin resolution
cp /projects/0/dfguu/data/hydroworld/PCRGLOBWB20/input5min/routing/lddsound_05min.map .
resample --clone lddsound_05min.map landmask_forcing.map landmask_forcing_05min.map
mapattr -c lddsound_05min.map landmask_forcing_05min.map
mapattr -p *.map
aguila landmask_forcing*
# identify all river basins
pcrcalc catchment_lddsound_05min.map = "catchment(lddsound_05min.map, pit(lddsound_05min.map))"
aguila catchment_lddsound_05min.map
# identify river basins with forcing - this will be the landmask for the extreme value analysis
pcrcalc areatotal_scalar_landmask_forcing_catchment_lddsound_05min.map = "areatotal(cover(scalar(landmask_forcing_05min.map), 0.0), catchment_lddsound_05min.map)"
# - only river basins with their 75% cells identified in the landmask of forcing
pcrcalc areatotal_scalar_all_catchment_lddsound_05min.map = "areatotal(cover(scalar(landmask_forcing_05min.map), 1.0), catchment_lddsound_05min.map)"
mapattr -p areatotal*
aguila areatotal*
pcrcalc landmask_extreme_value_analysis_05min.map = "if(areatotal_scalar_landmask_forcing_catchment_lddsound_05min.map gt (0.75 * areatotal_scalar_all_catchment_lddsound_05min.map), boolean(1.0))"
pcrcalc landmask_extreme_value_analysis_05min.map = "if(defined(lddsound_05min.map), landmask_extreme_value_analysis_05min.map) "
mapattr -p landmask_extreme_value_analysis_05min.map
aguila landmask_extreme_value_analysis_05min.map areatotal*