Commit 8d8dfe3d authored by Nick Bearson's avatar Nick Bearson
Browse files

added the ifflw command-line script for adding lw to iff

parent 35a1dba1
......@@ -10,6 +10,9 @@ import os
import numpy as np
from pyhdf.SD import *
import logging
log = logging.getLogger(__name__)
FLT_EPSILON_DEF = 1.19209290E-7
class DemLW(object):
......@@ -88,7 +91,7 @@ class DemLW(object):
# Retrieve the lw dataset from this subgrid
# print "Retrieving lw from %s" % (subgrid_path)
log.debug("Retrieving lw from %s" % (subgrid_path))
sd = SD(subgrid_path)
ds ="LandWater")
#!/usr/bin/env python
# encoding: utf-8
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:
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,
Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
from iocrutch import iocrutch
import argparse
import shutil
import numpy as np
import logging
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
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)
if args.output is not None:
shutil.copyfile(args.input, args.output)
o = args.output
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)
if __name__ == "__main__":
\ No newline at end of file
......@@ -2,11 +2,13 @@
from setuptools import setup, find_packages
description='Package for creating DEM Land/Water fields',
author='Nick Bearson',
entry_points={'console_scripts': ['ifflw=demlw.ifflw:main'],}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment