Page MenuHomec4science
Files F111317962

GVA_1301-13_LV95_shading_relief.sh
No OneTemporary
Actions

File Metadata

Created
Wed, Apr 30, 18:29

GVA_1301-13_LV95_shading_relief.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_GVA_50cm
ROOFTOP_MASK=SON_GVA_50cm_all
OUTPUT_BASE=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
######################################################
echo STARTING AT $(date)
TMP_LAYER_FILE=TMP_layers.txt
## 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_RAW=RAW_$OUTPUT_BASE\_$MONTH\_$HOUR\h
RELIEF_MASKED_NORM=$OUTPUT_BASE\_$MONTH\_$HOUR\h
EXPRESSION="$RELIEF_MASKED_NORM = float(max($RELIEF_RAW,0))/255.0*($ROOFTOP_MASK*0+1)"
## Execute
r.relief input=$ELEVATION output=$RELIEF_RAW altitude=$ALTITUDE azimuth=$AZIMUTH
r.mapcalc expression="$EXPRESSION" region=$REGION
## d.remove type=raster name=$HORIZON_RAW_OUT
echo Computed visible points and applied mask. Saved in layer $RELIEF_MASKED_NORM
echo $RELIEF_MASKED_NORM >> $TMP_LAYER_FILE
echo $MONTH $HOUR $RELIEF_MASKED_NORM >> $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