changed demlw aggregation to happen all in-code, copying what's done in PGE01

demlw/ifflw.py

@@ -6,82 +6,124 @@ Script for adding a land/sea mask to IFF
 Requires you to set the following env variables (or replace their defaults below):
 	DEMLW_DIR - the location of the SDPTK DEM files, available at:
 					ftp://edhs1.gsfc.nasa.gov/edhs/sdptk/DEMdata/
-	LW_CACHE_FILE - the location of lw-cache.npy; a cache file generated for this package which is used to speed up aggregation
 
 Author: Nick Bearson, nickb@ssec.wisc.edu
 Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
 """
 
-from iocrutch import iocrutch
-
 import argparse
 import shutil
-import numpy as np
 import logging
+
+from iocrutch import iocrutch
+import numpy as np
+
 log = logging.getLogger(__name__)
 
 import pkg_resources
 
 import os
+
 DEMLW_DIR = os.getenv("DEMLW_DIR", "anc/demlw/")
-LW_CACHE_FILE = os.getenv("LW_CACHE_FILE", os.path.dirname(__file__) + "/lw-cache.npy")
 
 from demlw import DemLW
 
+
 def aggregation(lsm):
-	"""
-	poor man's replication of PGE01/MOD_PR03/GEO_aggregate.c
-	while C6 looks at the 250m pixels inside a 1km footprint, we only
-	have the 1km pixels so we look at the ones surrounding it
-	"""
-	# FIXME: possible alternative - can we determine all LSM points that lie within a 1km footprint?
-	from scipy.ndimage.filters import generic_filter
-	xs, ys = lsm.shape
-	cache = np.load(LW_CACHE_FILE) # rather than compute all the combinations as they happen, we precomputed them and saved the lookup array
-	def sample(arr):
-		counts = tuple(np.bincount(arr.astype(np.uint8), minlength=8))
-		return cache[counts]
-	blurred = generic_filter(lsm, sample, 3)  # passes 3x3 box to sample for each point and assigns result
-	return blurred
+   """
+   poor man's replication of PGE01/MOD_PR03/GEO_aggregate.c
+   while C6 looks at the 250m pixels inside a 1km footprint, we only
+   have the 1km pixels so we look at the ones surrounding it
+   """
+   xs, ys = lsm.shape
+   newlsm = lsm.copy()
+   for x in range(xs):
+      for y in range(ys):
+         if x == 0 or x == xs - 1 or y == 0 or y == ys - 1: # can't do anything on the edges
+            continue
+         sample = lsm[x - 1:x + 2, y - 1:y + 2]  # sample a 3x3 box
+         newlsm[x][y] = aggregate_sample(sample)
+   return newlsm
+
+
+def aggregate_sample(sample):
+   SHALLOW_OCEAN, DRYLAND, COAST, SHALLOW_INLAND, EPHEMERAL, DEEP_INLAND, CONTINENTAL, DEEP_OCEAN = range(8)
+   counts = np.bincount(sample.flatten(), minlength=8)
+
+   land = counts[DRYLAND] + counts[COAST]
+
+   inland = counts[SHALLOW_INLAND] + counts[DEEP_INLAND]
+   ocean = counts[SHALLOW_OCEAN] + counts[CONTINENTAL] + counts[DEEP_OCEAN]
+   moderate = counts[DEEP_INLAND] + counts[CONTINENTAL]
+   water = inland + ocean + counts[EPHEMERAL]
+
+   lw_sum = land + water
+
+   if counts[EPHEMERAL] > (lw_sum / 2):
+      return EPHEMERAL
+   if (land + counts[EPHEMERAL]) > (lw_sum / 2):
+      if water > 0: # convert to coast if water > 0
+         return COAST
+      else:
+         return DRYLAND
+   else:
+      shallow = counts[EPHEMERAL] + counts[SHALLOW_INLAND] + counts[SHALLOW_OCEAN]
+      if shallow > (water / 2):
+         if inland > (water / 2):
+            return SHALLOW_INLAND
+         else:
+            return SHALLOW_OCEAN
+
+      if inland > (water / 2):
+         return DEEP_INLAND
+
+      shallow += moderate  # C:  /* Redefine "shallow" downwards */
+      if shallow > (water / 2):
+         return CONTINENTAL
+      else:
+         return DEEP_OCEAN
+
 
 def main():
-	parser = argparse.ArgumentParser(description="Utility for adding a land/sea mask to an Imager File Format (IFF) file.")
-	parser.add_argument('input', help='input file (full path)')
-	parser.add_argument('-a', '--aggregation', dest='agg', action="store_true", default=False, help='perform aggregation (MODIS C6-lite)')
-	parser.add_argument('-o', '--output', dest='output', help='output file (if you don\'t want to process in-place)')
-	parser.add_argument('-v', '--verbose', dest='verbosity', action="count", default=0,
-							  help='each occurrence increases verbosity through ERROR, WARNING, INFO, DEBUG')
-	args = parser.parse_args()
-
-	levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
-	verbosity = min(args.verbosity, len(levels)-1)
-	logging.basicConfig(level=levels[verbosity])
-
-	if args.output is not None:
-		shutil.copyfile(args.input, args.output)
-		o = args.output
-	else:
-		o = args.input
-
-
-	io = iocrutch.opensesame(o)
-	lats = io.get_variable_object("Latitude")[:]
-	lons = io.get_variable_object("Longitude")[:]
-
-	log.debug("Creating demlw object...")
-	demlw = DemLW(DEMLW_DIR)
-	log.debug("Performing demlw resampling...")
-	lw = demlw.latlon2demlw(lons, lats)
-
-	version = pkg_resources.require("demlw")[0].version # FIXME: why does this break when it's at the top?
-
-	if args.agg:
-		log.debug("Aggregating demlw...")
-		lw = aggregation(lw)
-		version = version + "+aggregation"
-	io.create_new_variable("LandSeaMask", lw, compression=0, dimensions=['line', 'element'])
-	io.set_global_attribute("with_demlw", version)
-	io.close()
+   parser = argparse.ArgumentParser(
+      description="Utility for adding a land/sea mask to an Imager File Format (IFF) file.")
+   parser.add_argument('input', help='input file (full path)')
+   parser.add_argument('-a', '--aggregation', dest='agg', action="store_true", default=False,
+                       help='perform aggregation (MODIS C6-lite)')
+   parser.add_argument('-o', '--output', dest='output', help='output file (if you don\'t want to process in-place)')
+   parser.add_argument('-v', '--verbose', dest='verbosity', action="count", default=0,
+                       help='each occurrence increases verbosity through ERROR, WARNING, INFO, DEBUG')
+   args = parser.parse_args()
+
+   levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+   verbosity = min(args.verbosity, len(levels) - 1)
+   logging.basicConfig(level=levels[verbosity])
+
+   if args.output is not None:
+      shutil.copyfile(args.input, args.output)
+      o = args.output
+   else:
+      o = args.input
+
+   io = iocrutch.opensesame(o)
+   lats = io.get_variable_object("Latitude")[:]
+   lons = io.get_variable_object("Longitude")[:]
+
+   log.debug("Creating demlw object...")
+   demlw = DemLW(DEMLW_DIR)
+   log.debug("Performing demlw resampling...")
+   lw = demlw.latlon2demlw(lons, lats)
+
+   version = pkg_resources.require("demlw")[0].version  # FIXME: why does this break when it's at the top?
+
+   if args.agg:
+      log.debug("Aggregating demlw...")
+      lw = aggregation(lw)
+      version = version + "+aggregation"
+   io.create_new_variable("LandSeaMask", lw, compression=0, dimensions=['line', 'element'])
+   io.set_global_attribute("with_demlw", version)
+   io.close()
+
 
 if __name__ == "__main__":
-	main()
\ No newline at end of file
+   main()
\ No newline at end of file

setup.py

@@ -2,7 +2,7 @@
 from setuptools import setup, find_packages
 
 setup(name='demlw',
-      version='1.0',
+      version='1.0.1',
       description='Package for creating DEM Land/Water fields',
       author='Nick Bearson',
       author_email='nickb@ssec.wisc.edu',