Page MenuHomec4science
Files F111435702

#GVA_1301-13_LV03_shading_combined.sh#
No OneTemporary
Actions

File Metadata

Created
Thu, May 1, 15:32

#GVA_1301-13_LV03_shading_combined.sh#
View Options

#!/bin/bash
################## USER SETTINGS ##############
WORKDIR=/scratch/walch/grassgis/workdir
MERGE_PYTHON_CODE=/home/walch/code/energy-potential/DOM/merge_shading_stats.py
# MAP_PATH_MAC=/Users/alinawalch/Documents/grassgis/LSN/shade_horizon
# MAP_PATH_CLUSTER=/scratch/walch/grassgis/LSN/shade_horizon
SOLPOS_BASE=$WORKDIR\/GVA_SOLPOS_RELIEF/GVA_SOLPOS
LOCATION=GVA_1301-13_LV03
echo TESTING FOR FULL YEAR
echo
# Layer names
ELEVATION=DOM_2m
ROOFTOP_MASK=SON_raster
OUTPUT_BASE=COMBI
HORIZON_LAYER=shade_horizon
RELIEF_LAYER=shade_relief
##---------- SUMMING STATISTICS -----------------##
## https://grass.osgeo.org/grass74/manuals/r.series.html
MIN_VISIBILITY=$1
VIS_PERC=$2
SUMMARY_STATS_BASE=$WORKDIR\/$LOCATION\_SHADE_$OUTPUT_BASE\_VIS_$VIS_PERC
# ###################################################
## ------------ INPUT/OUTPUT FILES -------------##
# Solar positions:
SOLPOS_FILE=$SOLPOS_BASE\.txt
LAYERFILE=$SOLPOS_BASE\_YEAR_layers.txt
echo Saving names of output layers in $LAYERFILE
# ############ Set UP GRASS ENVIRONMENT #############
mkdir $SUMMARY_STATS_BASE
## Set region
g.region raster=$ELEVATION
REGION=current
## Add mapsets of relief and horizon to current mapset
g.mapsets -p mapset=$HORIZON_LAYER\,$RELIEF_LAYER operation=add
######################################################
echo STARTING AT $(date)
TMP_LAYER_FILE=TMP_layers.txt
rm $TMP_LAYER_FILE ## remove file just in case it already exists
rm $LAYERFILE
## for-loop from list: loop over all time stamps and create boolean masks indicating which pixels are visible
while read MONTH HOUR AZIMUTH ALTITUDE; do
echo Month: $MONTH\; Hour: $HOUR\; Azimuth: $AZIMUTH\; Altitude: $ALTITUDE
## Set up variables and commands
RELIEF=RELIEF\_$MONTH\_$HOUR\h
HORIZON=HORIZON\_$MONTH\_$HOUR\h
COMBI_MAP=$OUTPUT_BASE\_$MONTH\_$HOUR\h
EXPRESSION="$COMBI_MAP = $RELIEF * $HORIZON"
## Execute
r.mapcalc expression="$EXPRESSION" region=$REGION
## d.remove type=raster name=$HORIZON_RAW_OUT
echo Computed combined hillshade with cast shadows. Saved in layer $COMBI_MAP
echo $COMBI_MAP >> $TMP_LAYER_FILE
echo $MONTH $HOUR $COMBI_MAP >> $LAYERFILE
echo FINISHED ITERATION $MONTH\-$HOUR\h AT $(date)
done < $SOLPOS_FILE
echo FINISHED COMPUTING AND MASKING REILEFS AT $(date)
echo
#########################################################
# Sum together all layers created in the previous loop
FILE_OF_SUMS=$OUTPUT_BASE\_MEAN
r.series -nz file=$TMP_LAYER_FILE output=$FILE_OF_SUMS method=average
rm $TMP_LAYER_FILE
# copy file to edit null values (set all values below visibility threshold to null)
TMP_FILE_OF_SUMS=TMP_$FILE_OF_SUMS
g.copy raster=$FILE_OF_SUMS\,$TMP_FILE_OF_SUMS
r.null $TMP_FILE_OF_SUMS setnull=0-$MIN_VISIBILITY
# compute summary statistics: null_cells are 'fully shaded', non_null_cells are partly shaded
r.univar -t $TMP_FILE_OF_SUMS zones=$ROOFTOP_MASK output=$SUMMARY_STATS_BASE\/summary_stats.csv --overwrite
# calculate a new mask of rooftops, with only not fully shaded cells
UNSHADED_MASK=$ROOFTOP_MASK\_vis_$VIS_PERC\_unshaded
EXPRESSION="$UNSHADED_MASK = int($ROOFTOP_MASK*($TMP_FILE_OF_SUMS*0+1))"
r.mapcalc expression="$EXPRESSION" region=$REGION --overwrite
echo FINISHED COMPUTING FULLY SHADED CELLS AT $(date)
echo
#########################################################
# compute statistics for all partly shaded roofs
## for-loop from list:
while read MONTH HOUR LAYER; do
echo Month: $MONTH\; Hour: $HOUR\; Layer: $LAYER\
r.univar -t map=$LAYER zones=$UNSHADED_MASK output=$SUMMARY_STATS_BASE\/stats_$MONTH\_$HOUR\h.csv --overwrite
done < $LAYERFILE
g.remove -f type=raster pattern=TMP_*
rm $LAYERFILE
echo FINISHED COMPUTING PARTLY SHADED CELLS AT $(date)
## RUN PYTHON SCRIPT merge_summary_files
source activate py2
python $MERGE_PYTHON_CODE $SUMMARY_STATS_BASE
source deactivate py2
echo FINISHED MERGING STATISTICS AT $(date)

Event Timeline

c4science · Help