In [1]:
import ephem
import json
from shapely.geometry import shape, mapping, Point
import fiona
import gdal
import ogr, osr
import struct
import math
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import datetime
import pytz
from altair import *
%matplotlib inline
plt.rcParams['figure.figsize'] = (18, 6)
In [2]:
with open ('inputs/site.geojson') as f:
js = json.load(f)
s = shape(js['features'][0]['geometry'])
s
Out[2]:
In [3]:
dem = gdal.Open('inputs/dem/filled.tif')
In [4]:
site = ogr.Geometry(ogr.wkbPoint) # create an ogr geom instead of shapely
site.AddPoint(s.x, s.y)
sr = dem.GetProjection()
destSR = osr.SpatialReference()
inSRS_converter = osr.SpatialReference()
inSRS_converter.ImportFromWkt(sr)
inSRS_proj4 = inSRS_converter.ExportToProj4()
destSR.ImportFromProj4(inSRS_proj4)
srcSR = osr.SpatialReference()
srcSR.ImportFromProj4("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
srTrans = osr.CoordinateTransformation(srcSR,destSR)
site_reproj = site
site.Transform(srTrans)
print(site_reproj.ExportToWkt())
In [5]:
gt=dem.GetGeoTransform()
rb=dem.GetRasterBand(1)
def get_elev_at_point(geotransform, rasterband, pointgeom):
mx,my=pointgeom.GetX(), pointgeom.GetY()
#Convert from map to pixel coordinates.
#Only works for geotransforms with no rotation.
px = int((mx - geotransform[0]) / geotransform[1]) #x pixel
py = int((my - geotransform[3]) / geotransform[5]) #y pixel
intval=rasterband.ReadAsArray(px,py,1,1)[0]
return intval[0]
site_elev = get_elev_at_point(gt, rb, site_reproj)
site_elev
Out[5]:
Did the viewshed analysis in QGIS using https://github.com/zoran-cuckovic/QGIS-visibility-analysis
Options:
Kept working in qgis for speed, will convert to this notebook eventually
In QGIS the horizon was:
Then:
degrees(atan((horizon_elev-site_elev)/distance)Bearing calculation:
degrees(atan(($x_at(0)-site_lon)/($y_at(0)-site_lat))) + (180 *((($y_at(0)-site_lat) < 0) + ((($x_at(0)-site_lon) < 0 AND ($y_at(0)-site_lat) >0)*2)))
In [6]:
df = pd.read_csv('inputs/viewshed/horizon_dist.csv')
df_skyview = df.set_index('bearing_from_site')
df_skyview.plot(y='inclination_from_site', style='.', figsize=(18,6), xlim=(0, 360))
Out[6]:
In [7]:
pd.set_option('display.float_format', lambda x: '%.3f' % x)
df_skyview.loc[df_skyview['inclination_from_site'].idxmax()]
Out[7]:
In [8]:
panel = ephem.Observer()
panel.lon = str(s.x) # this needs to be a string. derp
panel.lat = str(s.y)
panel.elevation = site_elev+2 # +2meters stays consistent with viewshed analysis
# panel.date = '2017/05/17 20:30:00' # utc time
panel
Out[8]:
In [9]:
from pytz import all_timezones, common_timezones
local = pytz.timezone('America/Anchorage')
def get_sunangles(time_utc, observer):
observer.date = time_utc
solar = ephem.Sun(observer)
a = {}
a['az'] = math.degrees(solar.az)
a['alt'] = math.degrees(solar.alt)
return a
get_sunangles(datetime.datetime.utcnow(), panel)
df_sun = pd.DataFrame(index=pd.date_range('2018-01-01 00:00', '2018-12-31 23:59', freq='5min', tz=local))
samples_per_hour=12 # adjust based on freq
df_sun['ts_utc'] = (df_sun.index).tz_convert(pytz.timezone('UTC'))
df_sun['az'] = df_sun.apply(lambda row: get_sunangles(row['ts_utc'], panel)['az'], axis=1)
df_sun['alt'] = df_sun.apply(lambda row: get_sunangles(row['ts_utc'], panel)['alt'], axis=1)
In [10]:
# df_sun.alt.plot(figsize=(18,6))
In [11]:
skyview_bearings = df_skyview.index.tolist()
azimuths = df_sun.az.tolist()
all_bearings = list(skyview_bearings)
all_bearings.extend(x for x in azimuths if x not in all_bearings)
df_allbearings = df_skyview.reindex(all_bearings)
df_allbearings.sort_index(inplace=True)
df_allbearings = df_allbearings.inclination_from_site.interpolate().ffill().bfill()
In [12]:
# df_allbearings.plot(kind='area', figsize=(18,6))
In [13]:
df_viz = df_sun.join(df_allbearings, on='az')
df_viz['visible'] = False
df_viz.loc[df_viz['alt'] > df_viz['inclination_from_site'], 'visible'] = True
In [14]:
series_dailysum = df_viz.groupby(pd.TimeGrouper(freq='D'))['visible'].sum()/samples_per_hour
series_dailysum.plot(figsize=(18,6), style='.')
Out[14]:
In [15]:
series_monthlysum = df_viz.groupby(pd.TimeGrouper(freq='M'))['visible'].sum()/samples_per_hour
series_monthlysum.plot(kind='bar', figsize=(18,6))
Out[15]:
In [16]:
series_dailysum.to_csv(path='dailysum.csv', index=True, date_format='%Y-%m-%d', float_format='%.2f')
series_monthlysum.to_csv(path='monthlysum.csv', index=True, date_format='%B', float_format='%.2f')
In [17]:
series_monthlysum
Out[17]:
In [18]:
min_alt = 15.0
df_viz['effective'] = False
df_viz.loc[(df_viz['alt'] > min_alt), 'effective'] = True
df_viz['above_min_alt'] = False
df_viz.loc[(df_viz['visible'] == True) & (df_viz['effective'] == True), 'above_min_alt'] = True
series_abovemin_dailysum = df_viz.groupby(pd.TimeGrouper(freq='D'))['effective'].sum()/samples_per_hour
series_effec_dailysum = df_viz.groupby(pd.TimeGrouper(freq='D'))['above_min_alt'].sum()/samples_per_hour
series_effec_monthlysum = df_viz.groupby(pd.TimeGrouper(freq='M'))['above_min_alt'].sum()/samples_per_hour
series_effec_monthlysum
Out[18]:
In [19]:
df_dailysums = pd.concat([series_dailysum, series_effec_dailysum, series_abovemin_dailysum], axis=1)
df_dailysums.columns = ['Sun visible above terrain', 'Sun visible and effective', 'Sun above 15 degrees']
In [20]:
fig, ax = plt.subplots(1, 1, figsize=(18,6))
df_dailysums.plot(ax=ax)
ax.set(xlabel='Date',
ylabel='Hours per day')
ax.legend(loc=2)
fig.autofmt_xdate()
In [21]:
fig, ax = plt.subplots(1, 1, figsize=(18,6))
df_dailysums['Sun above 15 degrees'].sub(df_dailysums['Sun visible and effective'], axis=0).plot(ax=ax)
ax.set(Title='Hours effective solar blocked by terrain')
Out[21]:
In [ ]: