Initializing, please wait...

demlw/__init__.py

+from demlw import DemLW

demlw/demlw.py

+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Author: Nick Bearson, nickb@ssec.wisc.edu
+Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+"""
+
+import os
+
+import numpy as np
+from pyhdf.SD import *
+
+FLT_EPSILON_DEF = 1.19209290E-7
+
+class DemLW(object):
+  """
+  """
+
+  def __init__(self, lw_dir):
+    """
+    lw_dir is where DEM HDF files reside
+    """
+    self.DEM_15_DICT = {
+      'horiGrids' : ['W180', 'W120', 'W60', 'E0', 'E60', 'E120'],
+      'vertGrids' : ['N90', 'N45', 'N0', 'S45'],
+      'demPath'   : lw_dir + "/dem15ARC_%s%s.hdf",
+    
+      'pixPerDegree' : 240,
+      'nVertPixInSubgrid' : 10800,
+      'nHoriPixInSubgrid' : 14400,
+    }
+    
+    
+    self.DEM_30_DICT = {
+      'horiGrids' : ['W180', 'W60', 'E60'],
+      'vertGrids' : ['N90', 'N0'],
+      'demPath'   : lw_dir + "/dem30ARC_%s%s.hdf",
+    
+      'pixPerDegree' : 120,
+      'nVertPixInSubgrid' : 10800,
+      'nHoriPixInSubgrid' : 14400,
+    }
+
+
+
+  def _map_subgrids(self, longitudes, latitudes, demdict):
+    lw = np.zeros(shape=latitudes.shape, dtype=np.uint8)
+  
+    pixPerDegree = demdict['pixPerDegree']
+    nVertPixInSubgrid = demdict['nVertPixInSubgrid']
+    nHoriPixInSubgrid = demdict['nHoriPixInSubgrid']
+    vertGrids = demdict['vertGrids']
+    horiGrids = demdict['horiGrids']
+    demPath = demdict['demPath']
+  
+    # Convert lat/lon to pixel space
+    lat_offset = ((( 90 - latitudes ) * pixPerDegree) + FLT_EPSILON_DEF).astype(int)
+    lon_offset = (((180 + longitudes) * pixPerDegree) + FLT_EPSILON_DEF).astype(int)
+  
+    # Figure out which subgrid each lat/lon resides in
+    vert_subgrid = (lat_offset / nVertPixInSubgrid).astype(int)
+    hori_subgrid = (lon_offset / nHoriPixInSubgrid).astype(int)
+  
+    # Get each pixel's index within its subgrid
+    subgrid_index_row = (lat_offset % nVertPixInSubgrid).astype(int)
+    subgrid_index_col = (lon_offset % nHoriPixInSubgrid).astype(int)
+  
+    # FIXME: LOOK INTO GDgetpixvalues function in the SDPTK
+  
+    # For each subgrid...
+    for v in range(len(vertGrids)):
+      for h in range(len(horiGrids)):
+        subgrid_h = horiGrids[h]
+        subgrid_v = vertGrids[v]
+        subgrid_path = demPath % (subgrid_h, subgrid_v)
+  
+        # Get a mask defining which points lie within this subgrid
+        check_grid_vert = vert_subgrid == v
+        check_grid_hori = hori_subgrid == h
+        check_grid = check_grid_vert & check_grid_hori
+  
+        # Retrieve just the pixels we want
+        wanted_rows = subgrid_index_row[check_grid]
+        wanted_cols = subgrid_index_col[check_grid]
+  
+        # If we don't have any wanted pixels in this subgrid, skip it!
+        if wanted_rows.size == 0 and wanted_cols.size == 0:
+          continue
+  
+        # Retrieve the lw dataset from this subgrid
+#        print "Retrieving lw from %s" % (subgrid_path)
+        sd = SD(subgrid_path)
+        ds = sd.select("LandWater")
+  
+        # FIXME: can we get around reading in the whole array by using fancy indexing with pyhdf?
+        lm = ds[:]
+        lm = lm.astype(np.uint8) # Have to cast here because 30-arc files are int8 while 15-arc are uint8
+        lw[check_grid] = lm[wanted_rows,wanted_cols]
+  
+    return lw
+  
+  
+  def latlon2demlw(self, longitudes, latitudes):
+  
+#    # First fill in with 30 arc sec data
+#    lw = self._map_subgrids(longitudes, latitudes, self.DEM_30_DICT)
+#  
+#    # Get a mask of the fill values
+#    fills = (lw == 221)
+#  
+#    # Now replace those fill values with values from 30 arc file
+#    lw[fills] = self._map_subgrids(longitudes[fills], latitudes[fills], self.DEM_15_DICT)
+
+    # FIXME: fills seem to have changed, just use ARC15 for now as speedup seems negligable
+    lw = self._map_subgrids(longitudes, latitudes, self.DEM_15_DICT)
+
+    return lw
+
+
+
+def test():
+  lons = np.array(np.arange(-180, 181,  .1))
+  lats = np.array(np.arange( 90,  -91, -.1))
+  longitude, latitude = np.meshgrid(lons, lats)
+  
+  print "Calling latlon2demlw on a %s array." % (str(longitude.shape))
+  demlw = DemLW("landwater")
+  lm = demlw.latlon2demlw(longitude, latitude)
+
+  import matplotlib.pyplot as plt
+  plot = plt.imshow(lm)
+  plt.savefig('lm.png')
+
+
+if __name__ == "__main__":
+  test()

setup.py

+
+from setuptools import setup, find_packages
+
+setup(name='demlw',
+      version='0.1',
+      description='Package for creating DEM Land/Water fields',
+      author='Nick Bearson',
+      author_email='nickb@ssec.wisc.edu',
+      packages=find_packages('.'),
+      install_requires=['pyhdf',
+                        'numpy'],
+      )