Page MenuHomec4science
Files F112474550

grass_shading_relief.sh
No OneTemporary
Actions

File Metadata

Created
Sat, May 10, 20:19

grass_shading_relief.sh
View Options

#!/bin/bash
################## USER SETTINGS ##############
WORKDIR=/Users/alinawalch/Documents/EPFL/data/grassgis_test
# 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/SOLPOS
SUMMARY_STATS_BASE=$WORKDIR\/'LSN_SHADE_RELIEF'
MONTH=01
echo TESTING FOR MONTH $MONTH
echo
# Layer names
ELEVATION=DOM_LSN
ROOFTOP_MASK=rooftop_select_LSN
OUTPUT_BASE=RELIEF
##---------- SUMMING STATISTICS -----------------##
## https://grass.osgeo.org/grass74/manuals/r.series.html
MIN_VISIBILITY=0.1
# ###################################################
## ------------ INPUT/OUTPUT FILES -------------##
# Solar positions:
SOLPOS_FILE=$SOLPOS_BASE\_$MONTH\.txt
LAYERFILE=$SOLPOS_BASE\_$MONTH\_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 --overwrite
## 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
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 --overwrite
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\_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

Event Timeline

c4science · Help