Homec4science
Diffusion solar_potential (aw_solar) Browse

solarPOT/grass_scripts/#GVA_1301-13_LV95_grass_shading_horizon.sh#
91fa34538c9daw_solar

Last Change
Raw File

#GVA_1301-13_LV95_grass_shading_horizon.sh#
View Options

#!/bin/bash
################## USER SETTINGS ##############
# WORKDIR=/Users/alinawalch/Documents/EPFL/data/grassgis_test
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_HORIZON/SOLPOS
SUMMARY_STATS_BASE=$WORKDIR\/'GVA_1301-13_50cm_SHADE_HORIZON_VIS_0'
echo TESTING FOR FULL YEAR
echo
# Layer names
ELEVATION=DOM_GVA_50cm
ROOFTOP_MASK=SON_GVA_50cm_all
OUTPUT_BASE=HORIZON
##---------- HORIZON SETTINGS -----------------##
# https://grass.osgeo.org/grass74/manuals/r.horizon.html
# Settings
MAXDISTANCE=5000
DISTANCE_FACTOR=1.0 # sampling distance step coefficient
REGION=current
##BUFFER=
##---------- SUMMING STATISTICS -----------------##
## https://grass.osgeo.org/grass74/manuals/r.series.html
VISIBILITY_PERC=0
# ###################################################
## ------------ INPUT/OUTPUT FILES -------------##
# Solar positions:
SOLPOS_FILE=$SOLPOS_BASE\.txt
LAYERFILE=$SOLPOS_BASE\_YEAR_layers.txt
echo Saving names of output layers in $LAYERFILE
# List of directions
DIRECTIONS_FILE=$SOLPOS_BASE\_direction.txt
# ############ Set UP GRASS ENVIRONMENT #############
mkdir $SUMMARY_STATS_BASE
## Set region
g.region raster=$ELEVATION
#################################################
echo STARTING AT $(date)
# file that is used to store only the names of the layers, without month and hour tags
# this file is used as input to r.series
# LOOP 1: loop over all distinct directions and compute horizons
while read DIRECTION; do
echo Direction: $DIRECTION
## Set up variables and commands
HORIZON_RAW=RAW_$OUTPUT_BASE
## Execute
## r.horizon -d $ELEVATION output=$HORIZON_RAW direction=$DIRECTION step=0.0 maxdistance=$MAXDISTANCE distance=$DISTANCE_FACTOR
echo FINISHED DIRECTION $DIRECTION AT $(date)
done < $DIRECTIONS_FILE
echo FINISHED COMPUTING HORIZONS AT $(date)
echo
######################################################
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 DIRECTION ALTITUDE; do
echo Month: $MONTH\; Hour: $HOUR\; Direction: $DIRECTION\; Altitude: $ALTITUDE
## Set up variables and commands
HORIZON_RAW=RAW_$OUTPUT_BASE
HORIZON_RAW_OUT=$HORIZON_RAW\_$DIRECTION\_0
HORIZON_MASKED_BOOL=$OUTPUT_BASE\_$MONTH\_$HOUR\h
EXPRESSION="$HORIZON_MASKED_BOOL = ($ALTITUDE>$HORIZON_RAW_OUT)*($ROOFTOP_MASK*0+1)"
## Execute
r.mapcalc expression="$EXPRESSION" region=$REGION
## d.remove type=raster name=$HORIZON_RAW_OUT
echo Computed visible points and applied mask. Saved in layer $HORIZON_MASKED_BOOL
echo $HORIZON_MASKED_BOOL >> $TMP_LAYER_FILE
echo $MONTH $HOUR $HORIZON_MASKED_BOOL >> $LAYERFILE
done < $SOLPOS_FILE
echo FINISHED MASKING HORIZONS 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-$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_$MIN_VISIBILITY\_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)
c4science · Help