Saturday, 25 December 2010

Daemonizing CherryPy

It is documented but without a proper example I found it very hard to find out how to daemonize a CherryPy server. After some trial and error it proved to be not so hard at all.

The CherryPy documentation is a little haphazard at times and although the Daemonizer plugin is documented I found it a bit difficult to understand without any examples. There is a bit more available now in the new documentation but that is quite hidden so it never hurts to show an example:

cherrypy.process.plugins.Daemonizer(cherrypy.engine).subscribe()
...
...
cherrypy.quickstart(Root(),config={
  '/':
  { 'log.access_file' : os.path.join(current_dir,"access.log"),
  'log.screen': False,
  'tools.sessions.on': True
  }})
Note that this only works on UNIX like systems (it certainly won't work on Windows XP). Also note the configuration parameters. Make sure sure you log your accesses explicitly otherwise you will have a hard time finding where the logging of your daemonized process went (hint: probably nowhere...)

Saturday, 18 December 2010

Calulating sunrise and sunset in Python

I expected to find dozens of readily available implementations of sunrise and sunset calculations in Python on the web but this turned out to be a disappointment. Therefore I resolved to write my own straightforward implementation.

The sunrise equation is pretty easy to find on Wikipedia but actual implementations are not, certainly not in Python 3.x. (If you are willing to stick to Python 2.x there is of course the excellent PyEphem package, see the end of this article) Fortunately NOAA provides annotated equations in the form of an OpenOffice spreadsheet. This allows for simple re-engineering in Python and gives us a way to verify the results of a new implementation against those in the spreadsheet.

If you save the code below as sunrise.py then calculating the time of sunrise today would be straightforward as shown in this example:

 import datetime
 import sunrise
 s = sun(lat=49,long=3)
 print('sunrise at ',s.sunrise(when=datetime.datetime.now())

The sun class also provides a sunset() method and solarnoon() method. All three methods take a when parameter that should be a datetime.datetime object. If this object contains timezone information or daylight saving time information, this information is used when calculating the times of sunrise, sunset and the solar noon.

Note that if no when parameter is given, a default datetime is used that is initialized with a LocalTimezone object from the timezone module. I have not provided that module here but you can implement one simple enough by copying the example in Python's documentation or you can comment out the import statement below and always supply a when parameter.

from math import cos,sin,acos,asin,tan
from math import degrees as deg, radians as rad
from datetime import date,datetime,time

# this module is not provided here. See text.
from timezone import LocalTimezone

class sun:
 """ 
 Calculate sunrise and sunset based on equations from NOAA
 http://www.srrb.noaa.gov/highlights/sunrise/calcdetails.html

 typical use, calculating the sunrise at the present day:
 
 import datetime
 import sunrise
 s = sun(lat=49,long=3)
 print('sunrise at ',s.sunrise(when=datetime.datetime.now())
 """
 def __init__(self,lat=52.37,long=4.90): # default Amsterdam
  self.lat=lat
  self.long=long
  
 def sunrise(self,when=None):
  """
  return the time of sunrise as a datetime.time object
  when is a datetime.datetime object. If none is given
  a local time zone is assumed (including daylight saving
  if present)
  """
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.sunrise_t)
  
 def sunset(self,when=None):
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.sunset_t)
  
 def solarnoon(self,when=None):
  if when is None : when = datetime.now(tz=LocalTimezone())
  self.__preptime(when)
  self.__calc()
  return sun.__timefromdecimalday(self.solarnoon_t)
  
 @staticmethod
 def __timefromdecimalday(day):
  """
  returns a datetime.time object.
  
  day is a decimal day between 0.0 and 1.0, e.g. noon = 0.5
  """
  hours  = 24.0*day
  h      = int(hours)
  minutes= (hours-h)*60
  m      = int(minutes)
  seconds= (minutes-m)*60
  s      = int(seconds)
  return time(hour=h,minute=m,second=s)

 def __preptime(self,when):
  """
  Extract information in a suitable format from when, 
  a datetime.datetime object.
  """
  # datetime days are numbered in the Gregorian calendar
  # while the calculations from NOAA are distibuted as
  # OpenOffice spreadsheets with days numbered from
  # 1/1/1900. The difference are those numbers taken for 
  # 18/12/2010
  self.day = when.toordinal()-(734124-40529)
  t=when.time()
  self.time= (t.hour + t.minute/60.0 + t.second/3600.0)/24.0
  
  self.timezone=0
  offset=when.utcoffset()
  if not offset is None:
   self.timezone=offset.seconds/3600.0
  
 def __calc(self):
  """
  Perform the actual calculations for sunrise, sunset and
  a number of related quantities.
  
  The results are stored in the instance variables
  sunrise_t, sunset_t and solarnoon_t
  """
  timezone = self.timezone # in hours, east is positive
  longitude= self.long     # in decimal degrees, east is positive
  latitude = self.lat      # in decimal degrees, north is positive

  time  = self.time # percentage past midnight, i.e. noon  is 0.5
  day      = self.day     # daynumber 1=1/1/1900
 
  Jday     =day+2415018.5+time-timezone/24 # Julian day
  Jcent    =(Jday-2451545)/36525    # Julian century

  Manom    = 357.52911+Jcent*(35999.05029-0.0001537*Jcent)
  Mlong    = 280.46646+Jcent*(36000.76983+Jcent*0.0003032)%360
  Eccent   = 0.016708634-Jcent*(0.000042037+0.0001537*Jcent)
  Mobliq   = 23+(26+((21.448-Jcent*(46.815+Jcent*(0.00059-Jcent*0.001813))))/60)/60
  obliq    = Mobliq+0.00256*cos(rad(125.04-1934.136*Jcent))
  vary     = tan(rad(obliq/2))*tan(rad(obliq/2))
  Seqcent  = sin(rad(Manom))*(1.914602-Jcent*(0.004817+0.000014*Jcent))+sin(rad(2*Manom))*(0.019993-0.000101*Jcent)+sin(rad(3*Manom))*0.000289
  Struelong= Mlong+Seqcent
  Sapplong = Struelong-0.00569-0.00478*sin(rad(125.04-1934.136*Jcent))
  declination = deg(asin(sin(rad(obliq))*sin(rad(Sapplong))))
  
  eqtime   = 4*deg(vary*sin(2*rad(Mlong))-2*Eccent*sin(rad(Manom))+4*Eccent*vary*sin(rad(Manom))*cos(2*rad(Mlong))-0.5*vary*vary*sin(4*rad(Mlong))-1.25*Eccent*Eccent*sin(2*rad(Manom)))

  hourangle= deg(acos(cos(rad(90.833))/(cos(rad(latitude))*cos(rad(declination)))-tan(rad(latitude))*tan(rad(declination))))

  self.solarnoon_t=(720-4*longitude-eqtime+timezone*60)/1440
  self.sunrise_t  =self.solarnoon_t-hourangle*4/1440
  self.sunset_t   =self.solarnoon_t+hourangle*4/1440

if __name__ == "__main__":
 s=sun(lat=52.37,long=4.90)
 print(datetime.today())
 print(s.sunrise(),s.solarnoon(),s.sunset())

For people willing to stick to Python 2.x there is a simple and good alternative in the form of the PyEphem package. It can do a lot more than just calculating sunsise. An example is shown below.

import ephem
o=ephem.Observer()
o.lat='49'
o.long='3'
s=ephem.Sun()
s.compute()
print ephem.localtime(o.next_rising(s))

Thursday, 16 December 2010

Python Geospatial Development

Every once in a while a book appears that immediately captures your attention. Python Geospatial Development by Erik Westra is such a book and the coming time I am sure I'll have an enjoyable time reading it. I'll post a full review when I am finished.

Done: Check this article.

Monday, 13 December 2010

Python thread safe cache class

Every so often the need arises to create a thread safe cache solution. This is my stab at a simple yet fully functional implementation that maintains the essential dictionary semantics, is thread safe and has a fixed, configurable size, for example in a multithreaded http server like CherryPy.

Although many dictionary operation like getting an item are reported to be atomic and therefore thread safe, this is actually an implementation specific feature of the widely used CPython implementation. And even so, adding keys or iterating over the keys in the dictionary might not be thread safe at all. We must therefore use some sort of locking mechanism to ensure no two threads try to modify the cache at the same time. (For more information check this discussion.)

The Cache class shown here features a configurable size and if the number of entries is too big it removes the oldest entry. We do not have to maintain a explicit usage administration for that because we make use of the properties of the OrderedDict class which remembers the order in which keys are inserted and sports a popitem() method that will remove the first (or last) item inserted.

from collections import OrderedDict
from threading import Lock

class Cache:
    def __init__(self,size=100):
        if int(size)<1 :
            raise AttributeError('size < 1 or not a number')
        self.size = size
        self.dict = OrderedDict()
        self.lock = Lock()

    def __getitem__(self,key):
        with self.lock:
            return self.dict[key]

    def __setitem__(self,key,value):
        with self.lock:
            while len(self.dict) >= self.size:
                self.dict.popitem(last=False)
            self.dict[key]=value

    def __delitem__(self,key):
        with self.lock:
            del self.dict[key]

Due to the functionality of the OrderedDict class we use, the implementation is very concise. The __init__() method merely checks whether the size attribute makes any sense and creates an instance of an OrderedDict and a Lock.

The with statements used in the remaining methods wait for the acquisition of the lock and guarantee that the lock is released even if an exception is raised. The __getitem__() method merely tries to retrieve a value by trying the key on the ordered dictionary after acquiring a lock.

The __setitem__() method removes as many items within its while loop to reduce the size to below the preset amount and then adds the new value. The popitem() method of an OrderedDict removes the least recently added key/value pair if it's last argument is set to False.

The __delitem__() also merely passes on the control to the underlying dictionary. Together these methods allow for any instance of our Cache class to be used like any other dictionary as the example code below illustrates:

>>> from cache import Cache
>>> c=Cache(size=3)
>>> c['key1']="one"
0
>>> c['key2']="two"
1
>>> c['key3']="three"
2
>>> c['key4']="four"
3
>>> c['key4']
'four'
>>> c['key1']
Traceback (most recent call last):
  File "", line 1, in 
  File "cache.py", line 13, in __getitem__
    return self.dict[key]
KeyError: 'key1'

Of course this doesn't show off the thread safety but it does show that the semantics are pretty much like that of a regular dictionary. If needed this class even be extended with suitable iterators/view like keys() and items() but for most caches this probably isn't necessary.