SF.net SVN: matplotlib:[7571] trunk/toolkits/basemap

[<js...@us...>]

Revision: 7571
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=7571&view=rev
Author:   jswhit
Date:     2009-08-26 01:12:32 +0000 (Wed, 26 Aug 2009)

Log Message:
-----------
added nightshade method

Modified Paths:
--------------
    trunk/toolkits/basemap/Changelog
    trunk/toolkits/basemap/examples/daynight.py
    trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
    trunk/toolkits/basemap/setup.py

Added Paths:
-----------
    trunk/toolkits/basemap/lib/mpl_toolkits/basemap/solar.py

Modified: trunk/toolkits/basemap/Changelog
===================================================================
--- trunk/toolkits/basemap/Changelog	2009-08-25 20:06:13 UTC (rev 7570)
+++ trunk/toolkits/basemap/Changelog	2009-08-26 01:12:32 UTC (rev 7571)
@@ -1,3 +1,6 @@
+version 0.99.5 (not yet released)
+           * added 'nightshade' method to shade night regions on a map.
+           * added lonmin, lonmax instance variables.
 version 0.99.4 (svn revision 7332)
            * ax.frame replaced with ax.spines to maintain compatibility 
              with matplotlib spines support.

Modified: trunk/toolkits/basemap/examples/daynight.py
===================================================================
--- trunk/toolkits/basemap/examples/daynight.py	2009-08-25 20:06:13 UTC (rev 7570)
+++ trunk/toolkits/basemap/examples/daynight.py	2009-08-26 01:12:32 UTC (rev 7571)
@@ -1,93 +1,16 @@
 import numpy as np
-from mpl_toolkits.basemap import Basemap, netcdftime
+from mpl_toolkits.basemap import Basemap
 import matplotlib.pyplot as plt
 from datetime import datetime
-from numpy import ma
 
 # example showing how to compute the day/night terminator and shade nightime
 # areas on a map.
-
-def epem(date):
-    """
-    input: date - datetime object (assumed UTC)
-    ouput: gha - Greenwich hour angle, the angle between the Greenwich
-           meridian and the meridian containing the subsolar point.
-           dec - solar declination.
-    """
-    dg2rad = np.pi/180.
-    rad2dg = 1./dg2rad
-    # compute julian day from UTC datetime object.
-    jday = netcdftime.JulianDayFromDate(date)
-    jd = np.floor(jday) # truncate to integer.
-    # utc hour.
-    ut = date.hour + date.minute/60. + date.second/3600.
-    # calculate number of centuries from J2000
-    t = (jd + (ut/24.) - 2451545.0) / 36525.
-    # mean longitude corrected for aberration
-    l = (280.460 + 36000.770 * t) % 360
-    # mean anomaly
-    g = 357.528 + 35999.050 * t
-    # ecliptic longitude
-    lm = l + 1.915 * np.sin(g*dg2rad) + 0.020 * np.sin(2*g*dg2rad)
-    # obliquity of the ecliptic
-    ep = 23.4393 - 0.01300 * t
-    # equation of time
-    eqtime = -1.915*np.sin(g*dg2rad) - 0.020*np.sin(2*g*dg2rad) \
-            + 2.466*np.sin(2*lm*dg2rad) - 0.053*np.sin(4*lm*dg2rad)
-    # Greenwich hour angle
-    gha = 15*ut - 180 + eqtime
-    # declination of sun
-    dec = np.arcsin(np.sin(ep*dg2rad) * np.sin(lm*dg2rad)) * rad2dg
-    return gha, dec
-
-def daynightgrid(date, nlons):
-    """
-    date is datetime object (assumed UTC).
-    nlons is # of longitudes used to compute terminator."""
-    dg2rad = np.pi/180.
-    lons = np.linspace(-180,180,nlons).astype(np.float32)
-    # compute greenwich hour angle and solar declination
-    # from datetime object (assumed UTC).
-    tau, dec = epem(date) 
-    # compute day/night terminator from hour angle, declination.
-    longitude = lons + tau
-    lats = np.arctan(-np.cos(longitude*dg2rad)/np.tan(dec*dg2rad))/dg2rad
-    # create day/night grid (1 for night, 0 for day)
-    nlats = ((nlons-1)/2)+1
-    lons2 = np.linspace(-180,180,nlons).astype(np.float32)
-    lats2 = np.linspace(-90,90,nlats).astype(np.float32)
-    lons2, lats2 = np.meshgrid(lons2,lats2)
-    lats = lats[np.newaxis,:]*np.ones((nlats,nlons),dtype=np.float32)
-    daynight = np.ones(lons2.shape, np.int8)
-    daynight = np.where(lats2>lats,0,daynight)
-    daynight = ma.array(daynight,mask=1-daynight) # mask day areas.
-    return lons2,lats2,daynight
-
-class Basemap2(Basemap):
-    def nightshade(self,date,color="k",nlons=1441,alpha=0.5,zorder=None):
-        # create grid of day=0, night=1
-        # 1441 means use a 0.25 degree grid.
-        lons,lats,daynight = daynightgrid(date,nlons)
-        x,y = self(lons,lats)
-        # contour the day-night grid, coloring the night area
-        # with the specified color and transparency.
-        CS = map.contourf(x,y,daynight,1,colors=[color],alpha=alpha)
-        # set zorder on ContourSet collections show night shading
-        # is on top.
-        for c in CS.collections:
-            if zorder is None:
-               c.set_zorder(c.get_zorder()+1)
-            else:
-               c.set_zorder(zorder)
-        return CS
-
-
 # miller projection 
-map = Basemap2(projection='mill',lon_0=0)
+map = Basemap(projection='mill',lon_0=180)
 # plot coastlines, draw label meridians and parallels.
 map.drawcoastlines()
 map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
-map.drawmeridians(np.arange(-180,180,60),labels=[0,0,0,1])
+map.drawmeridians(np.arange(map.lonmin,map.lonmax+30,60),labels=[0,0,0,1])
 # fill continents 'coral' (with zorder=0), color wet areas 'aqua'
 map.drawmapboundary(fill_color='aqua')
 map.fillcontinents(color='coral',lake_color='aqua')

Modified: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2009-08-25 20:06:13 UTC (rev 7570)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/__init__.py	2009-08-26 01:12:32 UTC (rev 7571)
@@ -53,7 +53,7 @@
 else:
     basemap_datadir = os.sep.join([os.path.dirname(__file__), 'data'])
 
-__version__ = '0.99.4'
+__version__ = '0.99.5'
 
 # supported map projections.
 _projnames = {'cyl'      : 'Cylindrical Equidistant',
@@ -3421,6 +3421,42 @@
             raise KeyError("barstyle must be 'simple' or 'fancy'")
         return rets
 
+    def nightshade(self,date,color="k",delta=0.25,alpha=0.5,ax=None,zorder=2):
+        """
+        Shade the regions of the map that are in darkness at the time
+        specifed by ``date``.  ``date`` is a datetime instance,
+        assumed to be UTC.
+
+        .. tabularcolumns:: |l|L|
+
+        ==============   ====================================================
+        Keywords         Description
+        ==============   ====================================================
+        color            color to shade night regions (default black).
+        delta            day/night terminator is computed with a
+                         a resolution of ``delta`` degrees (default 0.25).
+        alpha            alpha transparency for shading (default 0.5, so
+                         map background shows through).
+        zorder           zorder for shading (default 2).
+        ==============   ====================================================
+
+        Extra keyword ``ax`` can be used to override the default axis instance.
+
+        returns a matplotlib.contour.ContourSet instance.
+        """
+        from solar import daynight_grid
+        # create grid of day=0, night=1
+        lons,lats,daynight = daynight_grid(date,delta,self.lonmin,self.lonmax)
+        x,y = self(lons,lats)
+        # contour the day-night grid, coloring the night area
+        # with the specified color and transparency.
+        CS = self.contourf(x,y,daynight,1,colors=[color],alpha=alpha,ax=ax)
+        # set zorder on ContourSet collections show night shading
+        # is on top.
+        for c in CS.collections:
+            c.set_zorder(zorder)
+        return CS
+
     def _check_ax(self):
         """
         Returns the axis on which to draw.

Added: trunk/toolkits/basemap/lib/mpl_toolkits/basemap/solar.py
===================================================================
--- trunk/toolkits/basemap/lib/mpl_toolkits/basemap/solar.py	                        (rev 0)
+++ trunk/toolkits/basemap/lib/mpl_toolkits/basemap/solar.py	2009-08-26 01:12:32 UTC (rev 7571)
@@ -0,0 +1,66 @@
+# some simple functions to calculate solar position, day-night terminator
+import numpy as np
+from numpy import ma
+import netcdftime
+
+def epem(date):
+    """
+    input: date - datetime object (assumed UTC)
+    ouput: gha - Greenwich hour angle, the angle between the Greenwich
+           meridian and the meridian containing the subsolar point.
+           dec - solar declination.
+    """
+    dg2rad = np.pi/180.
+    rad2dg = 1./dg2rad
+    # compute julian day from UTC datetime object.
+    jday = netcdftime.JulianDayFromDate(date)
+    jd = np.floor(jday) # truncate to integer.
+    # utc hour.
+    ut = date.hour + date.minute/60. + date.second/3600.
+    # calculate number of centuries from J2000
+    t = (jd + (ut/24.) - 2451545.0) / 36525.
+    # mean longitude corrected for aberration
+    l = (280.460 + 36000.770 * t) % 360
+    # mean anomaly
+    g = 357.528 + 35999.050 * t
+    # ecliptic longitude
+    lm = l + 1.915 * np.sin(g*dg2rad) + 0.020 * np.sin(2*g*dg2rad)
+    # obliquity of the ecliptic
+    ep = 23.4393 - 0.01300 * t
+    # equation of time
+    eqtime = -1.915*np.sin(g*dg2rad) - 0.020*np.sin(2*g*dg2rad) \
+            + 2.466*np.sin(2*lm*dg2rad) - 0.053*np.sin(4*lm*dg2rad)
+    # Greenwich hour angle
+    gha = 15*ut - 180 + eqtime
+    # declination of sun
+    dec = np.arcsin(np.sin(ep*dg2rad) * np.sin(lm*dg2rad)) * rad2dg
+    return gha, dec
+
+def daynight_terminator(date, delta, lonmin, lonmax):
+    """
+    date is datetime object (assumed UTC).
+    nlons is # of longitudes used to compute terminator."""
+    dg2rad = np.pi/180.
+    lons = np.arange(lonmin,lonmax+0.5*delta,delta,dtype=np.float32)
+    # compute greenwich hour angle and solar declination
+    # from datetime object (assumed UTC).
+    tau, dec = epem(date) 
+    # compute day/night terminator from hour angle, declination.
+    longitude = lons + tau
+    lats = np.arctan(-np.cos(longitude*dg2rad)/np.tan(dec*dg2rad))/dg2rad
+    return lons, lats
+
+def daynight_grid(date, delta, lonmin, lonmax):
+    """
+    date is datetime object (assumed UTC).
+    delta is the grid interval (in degrees) used to compute terminator."""
+    lons, lats = daynight_terminator(date, delta, lonmin, lonmax)
+    # create day/night grid (1 for night, 0 for day)
+    lats2 = np.arange(-90,90+0.5*delta,delta,dtype=np.float32)
+    nlons = len(lons); nlats = len(lats2)
+    lons2, lats2 = np.meshgrid(lons,lats2)
+    lats = lats[np.newaxis,:]*np.ones((nlats,nlons),dtype=np.float32)
+    daynight = np.ones(lons2.shape, np.int8)
+    daynight = np.where(lats2>lats,0,daynight)
+    daynight = ma.array(daynight,mask=1-daynight) # mask day areas.
+    return lons2,lats2,daynight

Modified: trunk/toolkits/basemap/setup.py
===================================================================
--- trunk/toolkits/basemap/setup.py	2009-08-25 20:06:13 UTC (rev 7570)
+++ trunk/toolkits/basemap/setup.py	2009-08-26 01:12:32 UTC (rev 7571)
@@ -225,7 +225,7 @@
 package_data = {'mpl_toolkits.basemap':pyproj_datafiles+basemap_datafiles}
 setup(
   name              = "basemap",
-  version           = "0.99.4",
+  version           = "0.99.5",
   description       = "Plot data on map projections with matplotlib",
   long_description  = """
   An add-on toolkit for matplotlib that lets you plot data


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.