Astronomical libraries compatible with Ignition

DuctSoup · March 9, 2016, 11:25pm

Hi,

We use Ignition running on Linux to manage our buildings. I want to run a global gateway timer script to create and periodically update tags for things like:

[ul][li] local sunrise and sunset - power up or down the parking lot lights appropriately
[/li]
[li]sun position - to include with a transaction group for tracking the performance of solar PV and heating systems
[/li]
[li]theoretical solar radiation - for a local estimate of cloud cover to tip off our backup heat and lighting systems when the sun is taking a vacation day[/li][/ul]

I got really jazzed when I saw you could run 3rd party Python libraries. Turns out all the good libraries I’ve found so far are either written in or depend on C. I did get Pysolar to work but it’s incomplete, not being actively maintained and a little buggy.

Before I go reinvent the wheel or plug a Raspberry Pi in somewhere I thought I’d ask if anyone knows of a Python library for such things that’ll work with Ignition.

Thanks, Craig

Importing and using 3rd party Python libraries in Ignition
https://support.inductiveautomation.com/index.php?/Knowledgebase/Article/View/98/0/importing-and-using-3rd-party-python-libraries-in-ignition

Pysolar
http://pysolar.org/

Global Horizontal Irradiance Clear Sky Models: Implementation and Analysis
https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwiip_bN2bTLAhWrlIMKHT6sCaMQFgggMAA&url=http%3A%2F%2Fenergy.sandia.gov%2Fwp-content%2Fgallery%2Fuploads%2FSAND2012-2389_ClearSky_final.pdf&usg=AFQjCNHVHp2pbC-mQlPUL-iV-zIWpDrUig&sig2=iQdgRMlUcwQPagC7nsoRdw

diat150 · March 10, 2016, 12:53am

I dont know for sure but you may be able to grab some of this info from various web services across the internet. just a thought.

DuctSoup · March 10, 2016, 1:35am

That’s actually an excellent thought and exactly where we started a while back.

forecast-io-ignition
https://github.com/ductsoup/forecast-io-ignition

We do something similar with https://airnow.gov and a couple other JSON feeds. I didn’t bother publishing the rest to GitHub because they were just minor variations on the same idea.

The problem with the aggregate feeds is they’re dependent on our broadband connections, most don’t have timely data and none that I’ve seen have everything we want. That’s why we’re looking at doing the math on our Ignition gateways. Low overhead, higher reliability.

Kevin.Herron · March 10, 2016, 2:00am

You might look for libraries that do you what you need written in Java, which you can then integrate using the Ignition module SDK.

DuctSoup · March 10, 2016, 11:01pm

Yes, if necessary. A compatible Python library would be a more elegant and appropriate solution to this small problem.

DuctSoup · March 12, 2016, 2:45pm

Applied the 80/20 rule to a NOAA spreadsheet I found (that seems to have gone missing since the other day). It’s within a couple minutes for sunrise, noon and sunset and that’s close enough for what I’m doing.

Lots of room for improvement and not something I’d sign my name to at GitHub at this level. Thought I’d share here in case anybody else needed something basic.

Cut and paste the code into a shared gateway script library called something like noaa_solar
Edit for your lat, lon and tz and preferred root tag folder
Call it from a gateway timer script with “shared.noaa_solar.begin()” say about once a minute

[code]’’’
2016-03-12
http://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html
https://www.google.com/search?q=NOAA_Solar_Calculations_day&ie=utf-8&oe=utf-8#q=NOAA+Solar+Calculations+day
‘’’
import datetime, math

def begin():

def mkdir(par, dir):
	''' make a folder if necessary '''
	import system
	if (system.tag.exists(par + "/" + dir) == 0):
		system.tag.addTag(parentPath=par, name=dir, tagType="Folder")
	return

def mktag(par, tag, typ, val):
	''' make a tag if necessary then assign a value '''
	import system
	if (system.tag.exists(par + "/" + tag) == 0):
		system.tag.addTag(parentPath=par, name=tag, tagType="MEMORY", dataType=typ)
	system.tag.write(par + "/" + tag, val)
	return

def mkloc(par, lat, lon, tz):
	
	mktag(par, "Latitude", "Float8", lat)	
	mktag(par, "Longitude", "Float8", lon)
	mktag(par, "TimeZone", "Float8", tz)
	
	dt = datetime.datetime.now()
	#dt = datetime.datetime(2016, 3, 11, 0, 6, 0) #for testing
	mktag(par, "CurrentTime", "DateTime", dt)
	#print "date %s " % dt
	
	# E
	TimePastLocalMidnight = (dt - dt.replace(hour=0, minute=0, second=0, microsecond=0)).seconds / 60.0 / 1440
	Midnight = dt.replace(hour=0, minute=0, second=0, microsecond=0)
	mktag(par, "TimePastLocalMidnight", "Float8", TimePastLocalMidnight)
	#print "%s TimePastLocalMidnight"  % TimePastLocalMidnight
	
	# F 
	#ts = (dt - datetime.datetime(1970, 1, 1)).total_seconds()
	#revised for compatibility with Ignition
	dt = datetime.datetime(2016, 3, 11, 0, 6, 0) - datetime.datetime(1970, 1, 1)
	ts = dt.days * 86400.0 + dt.seconds
	JulianDay = ts / 86400.0 + 2440587.5 - tz / 24
	mktag(par, "JulianDay", "Float8", JulianDay)
	#print "%s JulianDay \n2457458.71" % JulianDay
	
	# G
	JulianCentury = (JulianDay - 2451545) / 36525
	mktag(par, "JulianCentury", "Float8", JulianCentury)
	#print "%s JulianCentury \n0.16190862" % JulianCentury
	
	# H
	GeomMeanLongSunDeg = 280.46646 + JulianCentury * (36000.76983 + JulianCentury * 0.0003032) % 360
	mktag(par, "GeomMeanLongSunDeg", "Float8", GeomMeanLongSunDeg)
	#print "%s GeomMeanLongSunDeg \n349.3015823" % GeomMeanLongSunDeg
	
	# I
	GeomMeanAnomSunDeg = 357.52911 + JulianCentury * (35999.05029 - 0.0001537 * JulianCentury)
	mktag(par, "GeomMeanAnomSunDeg", "Float8", GeomMeanAnomSunDeg)
	#print "%s GeomMeanAnomSunDeg \n6186.085812" % GeomMeanAnomSunDeg
	
	# J
	EccentEarthOrbit = 0.016708634 - JulianCentury * (0.000042037 + 0.0000001267 * JulianCentury)
	mktag(par, "EccentEarthOrbit", "Float8", EccentEarthOrbit)
	#print "%s EccentEarthOrbit \n0.016701825" % EccentEarthOrbit
	
	# K
	SunEqofCtr = math.sin(math.radians(GeomMeanAnomSunDeg)) * (1.914602 - JulianCentury * (0.004817 + 0.000014 * JulianCentury)) \
	   + math.sin(math.radians(2 * GeomMeanAnomSunDeg)) * (0.019993 - 0.000101 * JulianCentury) \
	   + math.sin(math.radians(3 * GeomMeanAnomSunDeg)) * 0.000289
	mktag(par, "SunEqofCtr", "Float8", SunEqofCtr)
	#print "%s SunEqofCtr \n1.764241916" % SunEqofCtr
	
	# L
	SunTrueLongDeg = GeomMeanLongSunDeg + SunEqofCtr
	mktag(par, "SunTrueLongDeg", "Float8", SunTrueLongDeg)
	#print "%s SunTrueLongDeg \n351.0658243" % SunTrueLongDeg
	
	# M
	SunTrueAnomDeg = GeomMeanAnomSunDeg + SunEqofCtr
	mktag(par, "SunTrueAnomDeg", "Float8", SunTrueAnomDeg)
	#print "%s SunTrueAnomDeg \n6187.850054" % SunTrueAnomDeg
	
	# N
	SunRadVectorAUs = (1.000001018 * (1 - EccentEarthOrbit * EccentEarthOrbit)) / (1 + EccentEarthOrbit * math.cos(math.radians(SunTrueAnomDeg)))
	mktag(par, "SunRadVectorAUs", "Float8", SunRadVectorAUs)
	#print "%s SunRadVectorAUs \n0.993466093" % SunRadVectorAUs
	
	# O
	SunAppLongDeg = SunTrueLongDeg - 0.00569 - 0.00478 * math.sin(math.radians(125.04 - 1934.136 * JulianCentury))
	mktag(par, "SunAppLongDeg", "Float8", SunAppLongDeg)
	#print "%s SunAppLongDeg \n351.0594597" % SunAppLongDeg
	
	# P
	MeanObliqEclipticDeg = 23 + (26 + ((21.448 - JulianCentury * (46.815 + JulianCentury * (0.00059 - JulianCentury * 0.001813))))/60)/60
	mktag(par, "MeanObliqEclipticDeg", "Float8", MeanObliqEclipticDeg)
	#print "%s MeanObliqEclipticDeg \n23.43718562" % MeanObliqEclipticDeg
	
	# Q
	ObliqCorrDeg = MeanObliqEclipticDeg + 0.00256 * math.cos(math.radians(125.04 - 1934.136 * JulianCentury))
	mktag(par, "ObliqCorrDeg", "Float8", ObliqCorrDeg)
	#print "%s ObliqCorrDeg \n23.43465125" % ObliqCorrDeg
	
	# R
	SunRtAscenDeg = math.degrees(math.atan2(math.cos(math.radians(SunAppLongDeg)),math.cos(math.radians(ObliqCorrDeg))*math.sin(math.radians(SunAppLongDeg)))) 
	mktag(par, "SunRtAscenDeg", "Float8", SunRtAscenDeg)
	#print "%s SunRtAscenDeg \n-8.213576963" % SunRtAscenDeg
	
	# S
	SunDeclinDeg = math.degrees(math.asin(math.sin(math.radians(ObliqCorrDeg))*math.sin(math.radians(SunAppLongDeg))))
	mktag(par, "SunDeclinDeg", "Float8", SunDeclinDeg)
	#print "%s SunDeclinDeg \n-3.543524949" % SunDeclinDeg
	
	# T
	vary = math.tan(math.radians(ObliqCorrDeg/2)) * math.tan(math.radians(ObliqCorrDeg/2))
	mktag(par, "vary", "Float8", vary)
	#print "%s vary \n0.043017009" % vary
	
	# U
	EqofTimeMin = 4 * math.degrees(vary * math.sin(2 * math.radians(GeomMeanLongSunDeg)) \
	   - 2 * EccentEarthOrbit * math.sin(math.radians(GeomMeanAnomSunDeg)) \
	   + 4 * EccentEarthOrbit * vary * math.sin(math.radians(GeomMeanAnomSunDeg)) * math.cos(2 * math.radians(GeomMeanLongSunDeg)) \
	   - 0.5 * vary * vary * math.sin(4 * math.radians(GeomMeanLongSunDeg)) \
	   - 1.25 * EccentEarthOrbit * EccentEarthOrbit * math.sin(2 * math.radians(GeomMeanAnomSunDeg)))
	mktag(par, "EqofTimeMin", "Float8", EqofTimeMin)
	#print "%s EqofTimeMin \n-9.949653921" % EqofTimeMin
	
	# V
	HASunriseDeg = math.degrees(math.acos(math.cos(math.radians(90.833)) / (math.cos(math.radians(lat)) * math.cos(math.radians(SunDeclinDeg))) \
	   - math.tan(math.radians(lat)) * math.tan(math.radians(SunDeclinDeg))))
	mktag(par, "HASunriseDeg", "Float8", HASunriseDeg)
	#print "%s HASunriseDeg \n88.11050567" % HASunriseDeg
	
	# W
	SolarNoonLST = (720 - 4 * lon - EqofTimeMin + tz * 60) / 1440
	#print Midnight + datetime.timedelta(seconds = 86400 * SolarNoonLST)
	mktag(par, "SolarNoonLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SolarNoonLST))
	#print "%s SolarNoonLST %s \n12:13:13" % (str(datetime.timedelta(seconds=SolarNoonLST * 86400)), SolarNoonLST)
	
	# X
	SunriseTimeLST = SolarNoonLST - HASunriseDeg * 4 / 1440
	mktag(par, "SunriseTimeLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SunriseTimeLST))
	#print "%s SunriseTimeLST %s \n6:20:47" % (str(datetime.timedelta(seconds=SunriseTimeLST * 86400)), SunriseTimeLST)
	
	# Y
	SunsetTimeLST = SolarNoonLST + HASunriseDeg * 4 / 1440
	mktag(par, "SunsetTimeLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SunsetTimeLST))
	#print "%s SunsetTimeLST %s \n18:05:40" % (str(datetime.timedelta(seconds=SunsetTimeLST * 86400)), SunsetTimeLST)
	
	# Z
	SunlightDurationMin = 8 * HASunriseDeg
	mktag(par, "SunlightDurationMin", "Float8", SunlightDurationMin)
	#print "%s SunlightDurationMin \n704.8840453" % SunlightDurationMin
	
	# AA
	TrueSolarTimeMin = (TimePastLocalMidnight * 1440 + EqofTimeMin + 4 * lon - 60 * tz) % 1440
	mktag(par, "TrueSolarTimeMin", "Float8", TrueSolarTimeMin)
	#print "%s TrueSolarTimeMin \n1432.778346" % TrueSolarTimeMin
	
	# AB
	#HourAngleDeg = IF(AA2/4<0,AA2/4+180,AA2/4-180)
	HourAngleDeg = TrueSolarTimeMin/4 + 180 if TrueSolarTimeMin/4 < 0 else TrueSolarTimeMin/4 - 180
	mktag(par, "HourAngleDeg", "Float8", HourAngleDeg)
	#print "%s HourAngleDeg \n178.1945865" % HourAngleDeg
	
	# AC
	SolarZenithAngleDeg = math.degrees(math.acos(math.sin(math.radians(lat)) * math.sin(math.radians(SunDeclinDeg)) + \
	   math.cos(math.radians(lat)) * math.cos(math.radians(SunDeclinDeg)) * math.cos(math.radians(HourAngleDeg))))
	mktag(par, "SolarZenithAngleDeg", "Float8", SolarZenithAngleDeg)
	#print "%s SolarZenithAngleDeg \n143.4909661" % SolarZenithAngleDeg
	
	# AD
	SolarElevationAngleDeg = 90 - SolarZenithAngleDeg
	mktag(par, "SolarElevationAngleDeg", "Float8", SolarElevationAngleDeg)
	#print "%s SolarElevationAngleDeg \n-53.49096606" % SolarElevationAngleDeg
	
	# AE
	if SolarElevationAngleDeg > 85:
	   ApproxAtmosphericRefractionDeg = 0,
	elif SolarElevationAngleDeg > 5:
	   ApproxAtmosphericRefractionDeg = 58.1 / math.tan(math.radians(SolarElevationAngleDeg)) - 0.07 / math.pow(math.tan(math.radians(SolarElevationAngleDeg)),3) + 0.000086 / math.pow(math.tan(math.radians(SolarElevationAngleDeg)),5)
	elif SolarElevationAngleDeg > -0.575:
	   ApproxAtmosphericRefractionDeg = 1735 + SolarElevationAngleDeg * (-518.2 + SolarElevationAngleDeg * (103.4 + SolarElevationAngleDeg * (-12.79 + SolarElevationAngleDeg * 0.711)))
	else:
	   ApproxAtmosphericRefractionDeg = -20.772 / math.tan(math.radians(SolarElevationAngleDeg))
	ApproxAtmosphericRefractionDeg /= 3600
	mktag(par, "ApproxAtmosphericRefractionDeg", "Float8", ApproxAtmosphericRefractionDeg)
	#print "%s ApproxAtmosphericRefractionDeg \n0.004270983" % ApproxAtmosphericRefractionDeg
	
	# AF
	SolarElevationCorrectedForATMRefractionDeg = SolarElevationAngleDeg + ApproxAtmosphericRefractionDeg
	mktag(par, "SolarElevationCorrectedForATMRefractionDeg", "Float8", SolarElevationCorrectedForATMRefractionDeg)
	#print "%s SolarElevationCorrectedForATMRefractionDeg \n-53.48669508" % SolarElevationCorrectedForATMRefractionDeg
	
	# AG
	if HourAngleDeg > 0:
	   SolarAzimuthAngleDegCWfromN = math.degrees(math.acos(((math.sin(math.radians(lat)) * math.cos(math.radians(SolarZenithAngleDeg))) \
	   - math.sin(math.radians(SunDeclinDeg))) / (math.cos(math.radians(lat)) * math.sin(math.radians(SolarZenithAngleDeg)))))+180 
	else:
	   SolarAzimuthAngleDegCWfromN = 540 - math.degrees(math.acos(((math.sin(math.radians(lat)) * math.cos(math.radians(SolarZenithAngleDeg))) \
	   - math.sin(math.radians(SunDeclinDeg))) / (math.cos(math.radians(lat)) * math.sin(math.radians(SolarZenithAngleDeg))))) 
	SolarAzimuthAngleDegCWfromN = SolarAzimuthAngleDegCWfromN % 360		
	mktag(par, "SolarAzimuthAngleDegCWfromN", "Float8", SolarAzimuthAngleDegCWfromN)
	#print "%s SolarAzimuthAngleDegCWfromN \n356.9703254" % SolarAzimuthAngleDegCWfromN
	#print "ok"
	
	return
	
root = "astro" # Ignition tag root folder

# check root exists  
mkdir("", root)

# create or update locations
mkloc(root + "/phl", 39.950, -75.167, -5)

return

[/code]