Skip to content
Browse files

initial commit with basic space gridding draft 1

  • Loading branch information...
1 parent d8efa22 commit 7de41a9cf8823c2bec4f06fff9a085919b59e494 Eva Schiffer committed Jan 21, 2014
Showing with 1,180 additions and 0 deletions.
  1. +5 −0 .gitignore
  2. +33 −0 setup.py
  3. 0 stg/__init__.py
  4. +35 −0 stg/constants.py
  5. +85 −0 stg/general_guidebook.py
  6. +131 −0 stg/io_manager.py
  7. +251 −0 stg/modis_guidebook.py
  8. +214 −0 stg/modis_io.py
  9. +20 −0 stg/rescale_filter.py
  10. +101 −0 stg/space_gridding.py
  11. +286 −0 stg/space_time_gridding.py
  12. +19 −0 stg/time_gridding.py
View
5 .gitignore
@@ -34,3 +34,8 @@ nosetests.xml
.mr.developer.cfg
.project
.pydevproject
+
+### Various OS Files ###
+.DS_Store
+._.DS_Store
+Thumbs.db
View
33 setup.py
@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+"""
+$Id$
+see http://peak.telecommunity.com/DevCenter/setuptools
+
+note: if you want to develop this code and run from code on the command line,
+please run the following line when you update to a new version of the code.
+python setup.py develop --install-dir=$HOME/Library/Python
+
+distribution:
+python setup.py develop --install-dir=$HOME/Library/Python
+python setup.py sdist
+python setup.py bdist_egg
+
+#(cd dist; rsync -Cuav * larch.ssec.wisc.edu:/var/apache/larch/htdocs/eggs/repos/uwglance)
+
+use:
+python setup.py install --install-dir=$HOME/Library/Python
+#easy_install -d $HOME/Library/Python -vi http://larch.ssec.wisc.edu/eggs/repos uwglance
+"""
+
+# changed to support egg distribution
+from setuptools import setup, find_packages
+
+setup( name="spacetimegrid",
+ version="0.1",
+ zip_safe = True,
+ entry_points = { 'console_scripts': [ 'stg = stg.space_time_gridding:main' ] },
+ packages = ['stg'], #find_packages('.'),
+ install_requires=[ 'numpy', 'scipy' ],
+ #package_data = {'': ['*.txt', '*.gif']}
+ )
+
View
0 stg/__init__.py
No changes.
View
35 stg/constants.py
@@ -0,0 +1,35 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Module to store all constants. Any constant needed in more than one
+component, or potentially more than one part should be defined here.
+
+Rules/Preferences:
+ - All values lowercase
+ - strings
+ - user-legible (assume that they may be printed in log messages)
+ - use == for comparison (not 'is' or 'not' or other)
+
+:author: Eva Schiffer (evas)
+:contact: eva.schiffer@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+
+Copyright (C) 2014 Space Science and Engineering Center (SSEC),
+ University of Wisconsin-Madison.
+"""
+__docformat__ = "restructuredtext en"
+
+
+
+SAT_AQUA = "aqua"
+SAT_TERRA = "terra"
+
+INST_MODIS = "modis"
+
+LON_KEY = "longitude"
+LAT_KEY = "latitude"
+DAY_MASK_KEY = "day_mask"
+NIGHT_MASK_KEY = "night_mask"
View
85 stg/general_guidebook.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Provide information about products the system expects, how they're stored
+in the files, and how to manage them during the process of regridding.
+
+:author: Eva Schiffer (evas)
+:contact: evas@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+:revision: $Id$
+"""
+__docformat__ = "restructuredtext en"
+
+from constants import *
+
+import sys
+import logging
+
+import stg.modis_guidebook as modis_guidebook
+
+LOG = logging.getLogger(__name__)
+
+def parse_datetime_from_filename (file_name_string) :
+ """parse the given file_name_string and create an appropriate datetime object
+ that represents the datetime indicated by the file name; if the file name does
+ not represent a pattern that is understood, None will be returned
+ """
+
+ datetime_to_return = None
+
+ # check to see the type of the file, then let the appropriate code parse it
+ if modis_guidebook.is_MODIS_file(file_name_string) :
+ datetime_to_return = modis_guidebook.parse_datetime_from_filename(file_name_string)
+
+ return datetime_to_return
+
+def get_satellite_from_filename (file_name_string) :
+ """given a file name, figure out which satellite it's from
+ if the file does not represent a known satellite name
+ configuration None will be returned
+ """
+
+ satellite_to_return = None
+ instrument_to_return = None
+
+ if modis_guidebook.is_MODIS_file(file_name_string) :
+ satellite_to_return, instrument_to_return = modis_guidebook.get_satellite_from_filename(file_name_string)
+
+ return satellite_to_return, instrument_to_return
+
+def get_variable_names (file_name_string, user_requested_names=[ ]) :
+ """get a list of variable names we expect to process from the file
+ """
+
+ var_names = set ( )
+
+ if modis_guidebook.is_MODIS_file(file_name_string) :
+ var_names.update(modis_guidebook.get_variable_names(user_requested_names))
+
+ return var_names
+
+def main():
+ import optparse
+ from pprint import pprint
+ usage = """
+%prog [options] filename1.hdf
+
+"""
+ parser = optparse.OptionParser(usage)
+ parser.add_option('-v', '--verbose', dest='verbosity', action="count", default=0,
+ help='each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG')
+ parser.add_option('-r', '--no-read', dest='read_hdf', action='store_false', default=True,
+ help="don't read or look for the hdf file, only analyze the filename")
+ (options, args) = parser.parse_args()
+
+ levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+ logging.basicConfig(level = levels[min(3, options.verbosity)])
+
+ LOG.info("Currently no command line tests are set up for this module.")
+
+if __name__ == '__main__':
+ sys.exit(main())
View
131 stg/io_manager.py
@@ -0,0 +1,131 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Handle parsing input from files.
+
+:author: Eva Schiffer (evas)
+:contact: evas@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+:revision: $Id$
+"""
+__docformat__ = "restructuredtext en"
+
+from constants import *
+
+import sys
+import logging
+import os
+
+import numpy
+
+import stg.general_guidebook as general_guidebook
+import stg.modis_guidebook as modis_guidebook
+import stg.modis_io as modis_io
+
+import keoni.fbf as fbf
+
+LOG = logging.getLogger(__name__)
+
+# these are suffixes used for temporary files
+DAY_TEMP_SUFFIX = "_daytemp"
+NIGHT_TEMP_SUFFIX = "_nighttemp"
+DAY_DENSITY_TEMP_SUFFIX = "_daydensitytemp"
+NIGHT_DENSITY_TEMP_SUFFIX = "_nightdensitytemp"
+EXPECTED_TEMP_SUFFIXES = [DAY_TEMP_SUFFIX, NIGHT_TEMP_SUFFIX, DAY_DENSITY_TEMP_SUFFIX, NIGHT_DENSITY_TEMP_SUFFIX]
+
+# these are suffixes used for the final, packed files
+DAY_SUFFIX = "_dayfinal"
+NIGHT_SUFFIX = "_nightfinal"
+EXPECTED_FINAL_SUFFIXES = [DAY_SUFFIX, NIGHT_SUFFIX]
+
+# all the suffixes we can produce
+ALL_EXPECTED_SUFFIXES = [DAY_TEMP_SUFFIX, NIGHT_TEMP_SUFFIX, DAY_DENSITY_TEMP_SUFFIX, NIGHT_DENSITY_TEMP_SUFFIX,
+ DAY_SUFFIX, NIGHT_SUFFIX]
+
+def open_file (file_path) :
+ """
+ given a file path that is a modis file, open it
+ """
+
+ file_object = None
+
+ if modis_guidebook.is_MODIS_file(file_path) :
+ file_object = modis_io.open_file(file_path)
+
+ return file_object
+
+def close_file (file_path, file_object) :
+ """
+ given a file object, close it
+ """
+
+ if modis_guidebook.is_MODIS_file(file_path) :
+ modis_io.close_file(file_object)
+
+def load_aux_data (file_path, minimum_scan_angle, file_object=None) :
+ """
+ load the auxillary data for a given file
+ """
+
+ temp_aux_data = None
+ if modis_guidebook.is_MODIS_file(file_path) :
+ file_object, temp_aux_data = modis_io.load_aux_data(file_path,
+ minimum_scan_angle,
+ file_object=file_object)
+
+ return file_object, temp_aux_data
+
+def load_variable_from_file (variable_name, file_path=None, file_object=None,
+ data_type_for_output=numpy.float32) :
+ """
+ load a given variable from a file path or file object
+ """
+
+ temp_data = None
+
+ if modis_guidebook.is_MODIS_file(file_path) :
+ file_object, temp_data = modis_io.load_variable_from_file (variable_name,
+ file_path=file_path,
+ file_object=file_object,
+ data_type_for_output=data_type_for_output)
+
+ return file_object, temp_data
+
+def save_data_to_file (stem_name, grid_shape, output_path, data_array, data_type, file_permissions="a") :
+ """
+ save numpy data to the appropriate file name given the information about the stem and path
+
+ by default this will append to the end of a file, you can instead pass in a different set of
+ file_permissions, using the standard permissions from python's open function
+ """
+
+ temp_file = fbf.filename(stem_name, data_type, shape=grid_shape)
+ temp_path = os.path.join(output_path, temp_file)
+ temp_file_obj = open(temp_path, 'a')
+ data_array.astype(data_type).tofile(temp_file_obj)
+ temp_file_obj.close()
+
+def main():
+ import optparse
+ from pprint import pprint
+ usage = """
+%prog [options] filename1.hdf
+
+"""
+ parser = optparse.OptionParser(usage)
+ parser.add_option('-v', '--verbose', dest='verbosity', action="count", default=0,
+ help='each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG')
+ parser.add_option('-r', '--no-read', dest='read_hdf', action='store_false', default=True,
+ help="don't read or look for the hdf file, only analyze the filename")
+ (options, args) = parser.parse_args()
+
+ levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+ logging.basicConfig(level = levels[min(3, options.verbosity)])
+
+ LOG.info("Currently no command line tests are set up for this module.")
+
+if __name__ == '__main__':
+ sys.exit(main())
View
251 stg/modis_guidebook.py
@@ -0,0 +1,251 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Provide information about products the system expects in MODIS files, how
+they're stored in the files, and how to manage them during the process of
+regridding.
+
+:author: Eva Schiffer (evas)
+:contact: evas@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+:revision: $Id$
+"""
+__docformat__ = "restructuredtext en"
+
+from constants import *
+
+import sys
+import os
+import logging
+import numpy
+from datetime import datetime
+
+LOG = logging.getLogger(__name__)
+
+# variable names expected in the files
+CLOUD_PHASE_NAME = 'Cloud_Phase_Infrared'
+CLOUD_TOP_TEMP_NAME = 'Cloud_Top_Temperature'
+CLOUD_TOP_PRESS_NAME = 'Cloud_Top_Pressure'
+CLOUD_EFF_EMISS_NAME = 'Cloud_Effective_Emissivity'
+CLOUD_TOP_PRESS_1KM_NAME = 'cloud_top_pressure_1km'
+CLOUD_TOP_TEMP_1KM_NAME = 'cloud_top_temperature_1km'
+CLOUD_EFF_RADIUS_16_NAME = 'Cloud_Effective_Radius_16'
+CLOUD_EFF_RADIUS_37_NAME = 'Cloud_Effective_Radius_37'
+CLOUD_OPTICAL_THICK_NAME = 'Cloud_Optical_Thickness'
+CLOUD_WATER_PATH_NAME = 'Cloud_Water_Path'
+QA_1KM_NAME = 'Quality_Assurance_1km'
+LATITUDE_NAME = 'Latitude'
+LONGITUDE_NAME = 'Longitude'
+SOLAR_ZENITH_NAME = 'Solar_Zenith'
+SENSOR_ZENITH_NAME = 'Sensor_Zenith'
+CLOUD_MULTI_LAYER_FLAG_NAME = 'Cloud_Multi_Layer_Flag'
+RADIANCE_VARIANCE_NAME = 'Radiance_Variance'
+BRIGHTNESS_TEMP_NAME = 'Brightness_Temperature'
+
+# TODO, sort out how to differentiate the Cloud Effective Radius
+# TODO, currently the fortran is testing:
+"""
+ if (versionnum<=4 & ~strncmp(filetype,'MAC',3))
+ varInFile = 'Effective_Particle_Radius';
+ else
+ varInFile = 'Cloud_Effective_Radius';
+ end
+"""
+CLOUD_EFF_RADIUS_NAME = ('Effective_Particle_Radius', 'Cloud_Effective_Radius')
+
+# important attribute names
+SCALE_ATTR_NAME = 'scale_factor'
+ADD_OFFSET_ATTR_NAME = 'add_offset'
+FILL_VALUE_ATTR_NAME = '_fillvalue'
+
+# a map of what the caller calls variables to variable names in files
+CALLER_VARIABLE_MAP = {
+ 'cloud phase': CLOUD_PHASE_NAME,
+ 'cloud_phase_infrared' : CLOUD_PHASE_NAME,
+ 'irphase': CLOUD_PHASE_NAME,
+
+ 'cloud top temperature': CLOUD_TOP_TEMP_NAME,
+ 'ctt5km': CLOUD_TOP_TEMP_NAME,
+ 'cloud_top_temperature': CLOUD_TOP_TEMP_NAME,
+
+ 'cloud top pressure': CLOUD_TOP_PRESS_NAME,
+ 'ctp5km': CLOUD_TOP_PRESS_NAME,
+ 'pressure': CLOUD_TOP_PRESS_NAME,
+
+ 'effective cloud emissivity': CLOUD_EFF_EMISS_NAME,
+
+ """ TODO, this won't work for cloud effective radius
+ 're': CLOUD_EFF_RADIUS_NAME,
+ 'effective_particle_radius': CLOUD_EFF_RADIUS_NAME,
+ 'effective radius': CLOUD_EFF_RADIUS_NAME,
+ 'cloud_effective_radius': CLOUD_EFF_RADIUS_NAME,
+ """
+
+ 'ctp1km': CLOUD_TOP_PRESS_1KM_NAME,
+
+ 'ctt1km': CLOUD_TOP_TEMP_1KM_NAME,
+
+ 're16': CLOUD_EFF_RADIUS_16_NAME,
+ 'cloud_effective_radius_16': CLOUD_EFF_RADIUS_16_NAME,
+
+ 're37': CLOUD_EFF_RADIUS_37_NAME,
+ 'cloud_effective_radius_37': CLOUD_EFF_RADIUS_37_NAME,
+
+ 'tau': CLOUD_OPTICAL_THICK_NAME,
+ 'cloud_optical_thickness': CLOUD_OPTICAL_THICK_NAME,
+ 'optical thickness': CLOUD_OPTICAL_THICK_NAME,
+
+ 'cwp': CLOUD_WATER_PATH_NAME,
+ 'cloud_water_path': CLOUD_WATER_PATH_NAME,
+ 'cloud water path': CLOUD_WATER_PATH_NAME,
+
+ 'retphase': QA_1KM_NAME,
+ 'retrieval phase': QA_1KM_NAME,
+ 'retqflag': QA_1KM_NAME,
+ 'retrieval_quality_flag': QA_1KM_NAME,
+
+ 'latitude': LATITUDE_NAME,
+ 'lat': LATITUDE_NAME,
+
+ 'longitude': LONGITUDE_NAME,
+ 'lon': LONGITUDE_NAME,
+
+ 'solar zenith': SOLAR_ZENITH_NAME,
+ 'sunzen': SOLAR_ZENITH_NAME,
+
+ 'sensor zenith': SENSOR_ZENITH_NAME,
+ 'satzen': SENSOR_ZENITH_NAME,
+ 'viewing zenith': SENSOR_ZENITH_NAME,
+
+ 'multilayer': CLOUD_MULTI_LAYER_FLAG_NAME,
+ 'cloud_multi_layer_flag': CLOUD_MULTI_LAYER_FLAG_NAME,
+ 'overlap': CLOUD_MULTI_LAYER_FLAG_NAME,
+
+ 'variance': RADIANCE_VARIANCE_NAME,
+ 'radiance_variance': RADIANCE_VARIANCE_NAME,
+ 'radiance_variability': RADIANCE_VARIANCE_NAME,
+
+ 'brightness_temperature': BRIGHTNESS_TEMP_NAME,
+ 'bt': BRIGHTNESS_TEMP_NAME,
+ }
+
+DATA_TYPE_TO_USE = { # TODO, eventually differentiate the data type by variable
+ CLOUD_PHASE_NAME: numpy.float32,
+ CLOUD_TOP_TEMP_NAME: numpy.float32,
+ CLOUD_TOP_PRESS_NAME: numpy.float32,
+ CLOUD_EFF_EMISS_NAME: numpy.float32,
+ CLOUD_TOP_PRESS_1KM_NAME: numpy.float32,
+ CLOUD_TOP_TEMP_1KM_NAME: numpy.float32,
+ CLOUD_EFF_RADIUS_16_NAME: numpy.float32,
+ CLOUD_EFF_RADIUS_37_NAME: numpy.float32,
+ CLOUD_OPTICAL_THICK_NAME: numpy.float32,
+ CLOUD_WATER_PATH_NAME: numpy.float32,
+ QA_1KM_NAME: numpy.float32,
+ LATITUDE_NAME: numpy.float32,
+ LONGITUDE_NAME: numpy.float32,
+ SOLAR_ZENITH_NAME: numpy.float32,
+ SENSOR_ZENITH_NAME: numpy.float32,
+ CLOUD_MULTI_LAYER_FLAG_NAME: numpy.float32,
+ RADIANCE_VARIANCE_NAME: numpy.float32,
+ BRIGHTNESS_TEMP_NAME: numpy.float32,
+ }
+
+# a list of the default variables expected in a file, used when no variables are selected by the caller
+EXPECTED_VARIABLES_IN_FILE = set([CLOUD_TOP_PRESS_NAME]) # TODO, this is currently set for our minimal testing
+
+# TODO, move this up to the general_guidebook
+def _clean_off_path_if_needed(file_name_string) :
+ """
+ remove the path from the file if nessicary
+ """
+
+ return os.path.basename(file_name_string)
+
+def is_MODIS_file (file_name_string) :
+ """determine if a file name is the right pattern to represent a MODIS file
+ if the file_name_string matches how we expect MODIS files to look return
+ TRUE else will return FALSE
+ """
+
+ temp_name_string = _clean_off_path_if_needed(file_name_string)
+
+ return (temp_name_string.startswith('MYD') or temp_name_string.startswith('MOD')) and temp_name_string.endswith('hdf')
+
+def parse_datetime_from_filename (file_name_string) :
+ """parse the given file_name_string and create an appropriate datetime object
+ that represents the datetime indicated by the file name; if the file name does
+ not represent a pattern that is understood, None will be returned
+ """
+
+ temp_name_string = _clean_off_path_if_needed(file_name_string)
+
+ datetime_to_return = None
+
+ # there are at least two file name formats to parse here
+ if temp_name_string.startswith('MYD') or temp_name_string.startswith('MOD') :
+ temp = temp_name_string.split('.')
+ datetime_to_return = datetime.strptime(temp[1] + temp[2], "A%Y%j%H%M")
+ # I confirmed with Nick that this is the correct date format
+
+ return datetime_to_return
+
+def get_satellite_from_filename (data_file_name_string) :
+ """given a file name, figure out which satellite it's from
+ if the file does not represent a known satellite name
+ configuration None will be returned
+ """
+
+ temp_name_string = _clean_off_path_if_needed(data_file_name_string)
+
+ satellite_to_return = None
+ instrument_to_return = None
+
+ if temp_name_string.startswith("MYD") :
+ satellite_to_return = SAT_AQUA
+ instrument_to_return = INST_MODIS
+ elif temp_name_string.startswith("MOD") :
+ satellite_to_return = SAT_TERRA
+ instrument_to_return = INST_MODIS
+
+ return satellite_to_return, instrument_to_return
+
+def get_variable_names (user_requested_names) :
+ """get a list of variable names we expect to process from the file
+ """
+
+ var_names = set( )
+
+ if len(user_requested_names) <= 0 :
+ var_names.update(EXPECTED_VARIABLES_IN_FILE)
+ else :
+
+ for user_name in user_requested_names :
+ if user_name in CALLER_VARIABLE_MAP.keys() :
+ var_names.update(set([CALLER_VARIABLE_MAP[user_name]]))
+
+ return var_names
+
+def main():
+ import optparse
+ from pprint import pprint
+ usage = """
+%prog [options] filename1.hdf
+
+"""
+ parser = optparse.OptionParser(usage)
+ parser.add_option('-v', '--verbose', dest='verbosity', action="count", default=0,
+ help='each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG')
+ parser.add_option('-r', '--no-read', dest='read_hdf', action='store_false', default=True,
+ help="don't read or look for the hdf file, only analyze the filename")
+ (options, args) = parser.parse_args()
+
+ levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+ logging.basicConfig(level = levels[min(3, options.verbosity)])
+
+ LOG.info("Currently no command line tests are set up for this module.")
+
+if __name__ == '__main__':
+ sys.exit(main())
View
214 stg/modis_io.py
@@ -0,0 +1,214 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+Handle parsing input from Modis files.
+
+:author: Eva Schiffer (evas)
+:contact: evas@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+:revision: $Id$
+"""
+__docformat__ = "restructuredtext en"
+
+from constants import *
+
+import sys
+import logging
+
+import numpy
+from pyhdf.SD import SD,SDC, SDS, HDF4Error
+
+from stg.constants import *
+import stg.modis_guidebook as modis_guidebook
+
+LOG = logging.getLogger(__name__)
+
+# the line between day and night for our day/night masks (in solar zenith angle degrees)
+DAY_NIGHT_LINE_DEGREES = 84.0
+
+def open_file (file_path) :
+ """
+ given a file path that is a modis file, open it
+ """
+
+ file_object = SD(file_path, SDC.READ)
+
+ return file_object
+
+def close_file (file_object) :
+ """
+ given a file object, close it
+ """
+
+ file_object.end()
+
+def load_aux_data (file_path, minimum_scan_angle, file_object=None) :
+ """
+ load the auxillary data and process the appropriate masks from it
+ """
+
+ # make our return structure
+ aux_data_sets = { }
+
+ # load the longitude and latitude
+ file_object, aux_data_sets[LON_KEY] = load_variable_from_file (modis_guidebook.LONGITUDE_NAME,
+ file_path=file_path, file_object=file_object)
+ file_object, aux_data_sets[LAT_KEY] = load_variable_from_file (modis_guidebook.LATITUDE_NAME,
+ file_path=file_path, file_object=file_object)
+
+ # load the angles to make masks
+ file_object, solar_zenith_data_temp = load_variable_from_file (modis_guidebook.SOLAR_ZENITH_NAME,
+ file_path=file_path, file_object=file_object)
+ file_object, sat_zenith_data_temp = load_variable_from_file (modis_guidebook.SENSOR_ZENITH_NAME,
+ file_path=file_path, file_object=file_object)
+
+ # transform the satellite zenith to scan angle
+ scan_angle_data_temp = satellite_zenith_angle_to_scan_angle(sat_zenith_data_temp)
+
+ # build the day and night masks
+ ok_scan_angle = scan_angle_data_temp <= minimum_scan_angle
+ aux_data_sets[DAY_MASK_KEY] = (solar_zenith_data_temp < DAY_NIGHT_LINE_DEGREES) & ok_scan_angle
+ aux_data_sets[NIGHT_MASK_KEY] = (solar_zenith_data_temp >= DAY_NIGHT_LINE_DEGREES) & ok_scan_angle
+
+ return file_object, aux_data_sets
+
+# FUTURE, the data type needs to be handled differently
+def load_variable_from_file (variable_name, file_path=None, file_object=None,
+ fill_value_name=modis_guidebook.FILL_VALUE_ATTR_NAME,
+ scale_name=modis_guidebook.SCALE_ATTR_NAME,
+ offset_name=modis_guidebook.ADD_OFFSET_ATTR_NAME,
+ data_type_for_output=numpy.float32) :
+ """
+ load a given variable from a file path or file object
+ """
+
+ if file_path is None and file_object is None :
+ raise ValueError("File path or file object must be given to load file.")
+ if file_object is None :
+ file_object = open_file(file_path)
+
+ variable_names = file_object.datasets().keys()
+ if variable_name not in variable_names :
+ raise ValueError("Variable " + str(variable_name) +
+ " is not present in file " + str(file_path) + " .")
+
+ # defaults
+ scale_factor = 1.0
+ add_offset = 0.0
+ data_type = None
+ scaling_method = None
+
+ # get the variable object and use it to
+ # get our raw data and scaling info
+ variable_object = file_object.select(variable_name)
+ raw_data_copy = variable_object[:]
+ raw_data_copy = raw_data_copy.astype(data_type_for_output) if data_type_for_output is not None else raw_data_copy
+ temp_attrs = variable_object.attributes()
+ try :
+ scale_factor, scale_factor_error, add_offset, add_offset_error, data_type = SDS.getcal(variable_object)
+ except HDF4Error:
+ # load just the scale factor and add offset information by hand
+ if offset_name in temp_attrs.keys() :
+ add_offset = temp_attrs[offset_name]
+ data_type = numpy.dtype(type(add_offset))
+ if scale_name in temp_attrs.keys() :
+ scale_factor = temp_attrs[scale_name]
+ data_type = numpy.dtype(type(scale_factor))
+
+ # get the fill value
+ fill_value = temp_attrs[fill_value_name] if fill_value_name in temp_attrs.keys() else numpy.nan
+ # change the fill value to numpy.nan
+ fill_mask = raw_data_copy == fill_value
+ if fill_value is not numpy.nan :
+ raw_data_copy[fill_mask] = numpy.nan
+ fill_value = numpy.nan
+
+ # we got all the info we need about that file
+ SDS.endaccess(variable_object)
+
+ # don't do lots of work if we don't need to scale things
+ if (scale_factor == 1.0) and (add_offset == 0.0) :
+ return file_object, raw_data_copy
+
+ # if we don't have a data type something strange has gone wrong
+ assert(data_type is not None)
+
+ # create the scaled version of the data
+ scaled_data_copy = raw_data_copy.copy()
+ scaled_data_copy = unscale_data(scaled_data_copy, fill_mask=fill_mask,
+ scale_factor=scale_factor, offset=add_offset)
+
+ return file_object, scaled_data_copy
+
+def unscale_data (data, fill_mask=None, scale_factor=None, offset=None) :
+ """
+ unscale the given data
+
+ data is modified in place and fill values will not be changed
+ if a scale factor or offset is given as None (or not given) it will not be applied
+
+ the general formula for unscaling data is:
+
+ final_data = scale_factor * (input_data - offset)
+
+ """
+
+ to_return = data
+
+ not_fill_mask = ~fill_mask
+
+ # if we have an offset use it to offset the data
+ if (offset is not None) and (offset is not 0.0) :
+ to_return[not_fill_mask] -= offset
+
+ # if we found a scale use it to scale the data
+ if (scale_factor is not None) and (scale_factor is not 1.0) :
+ to_return[not_fill_mask] *= scale_factor
+
+ return to_return
+
+# FUTURE, will this be used for other satellites? should it move up to the io_manager?
+def satellite_zenith_angle_to_scan_angle (sat_zenith_data) :
+ """
+ given a set of satellite zenith angles, calculate the equivalent scan angles
+
+ Note: This comes directly from Nadia's satz2scang function.
+ """
+
+ # some constants
+ re = 6371.03;
+ altkm = 825;
+ fac = re / (re + altkm);
+ dtr = 0.01745329;
+
+ # do the angle calculations
+ arg_data = sat_zenith_data * dtr
+ ang_data = numpy.sin(numpy.sin(arg_data) * fac)
+ scan_angle_data = ang_data / dtr
+
+ return scan_angle_data
+
+def main():
+ import optparse
+ from pprint import pprint
+ usage = """
+%prog [options] filename1.hdf
+
+"""
+ parser = optparse.OptionParser(usage)
+ parser.add_option('-v', '--verbose', dest='verbosity', action="count", default=0,
+ help='each occurrence increases verbosity 1 level through ERROR-WARNING-INFO-DEBUG')
+ parser.add_option('-r', '--no-read', dest='read_hdf', action='store_false', default=True,
+ help="don't read or look for the hdf file, only analyze the filename")
+ (options, args) = parser.parse_args()
+
+ levels = [logging.ERROR, logging.WARN, logging.INFO, logging.DEBUG]
+ logging.basicConfig(level = levels[min(3, options.verbosity)])
+
+ LOG.info("Currently no command line tests are set up for this module.")
+
+if __name__ == '__main__':
+ sys.exit(main())
View
20 stg/rescale_filter.py
@@ -0,0 +1,20 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+This module holds support functions that are used to rescale, filter, or
+otherwise transform the data before any gridding is done.
+
+:author: Eva Schiffer (evas)
+:contact: eva.schiffer@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+
+Copyright (C) 2014 Space Science and Engineering Center (SSEC),
+ University of Wisconsin-Madison.
+"""
+__docformat__ = "restructuredtext en"
+
+
+
View
101 stg/space_gridding.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+This module handles the basic space gridding algorithm.
+
+:author: Eva Schiffer (evas)
+:contact: eva.schiffer@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+
+Copyright (C) 2014 Space Science and Engineering Center (SSEC),
+ University of Wisconsin-Madison.
+"""
+__docformat__ = "restructuredtext en"
+
+import numpy
+import logging
+
+from stg.constants import *
+
+LOG = logging.getLogger(__name__)
+
+def calculate_index_from_nav_data (aux_data, grid_degrees) :
+ """
+ given the aux data, use the navigation and masks to calculate
+ where the elements will be space gridded
+ """
+
+ # figure out where the day/night indexes will fall
+ day_lon_index = numpy.round(aux_data[LON_KEY][aux_data[DAY_MASK_KEY]] / grid_degrees)
+ day_lat_index = numpy.round(aux_data[LAT_KEY][aux_data[DAY_MASK_KEY]] / grid_degrees)
+ night_lon_index = numpy.round(aux_data[LON_KEY][aux_data[NIGHT_MASK_KEY]] / grid_degrees)
+ night_lat_index = numpy.round(aux_data[LAT_KEY][aux_data[NIGHT_MASK_KEY]] / grid_degrees)
+
+ return day_lon_index, day_lat_index, night_lon_index, night_lat_index
+
+def space_grid_data (grid_lon_size, grid_lat_size, data, lon_indexes, lat_indexes ) :
+ """
+ given lon/lat indexes, data, and the grid size, sort the data into a space grid
+
+ returns the filled space grid (empty space is NaN values), a density map of where the data is, and the size of the deepest bucket
+ """
+
+ space_grid_shape = (grid_lon_size, grid_lat_size) # TODO, is this the correct order?
+
+ # create the density map and figure out how dense the data will be
+ # FUTURE, I do not like this looping solution, figure out how to do this in native numpy ops
+ density_map = numpy.zeros(space_grid_shape)
+ for index in range(len(data)) :
+ if not numpy.isnan(data[index]) :
+ density_map[lon_indexes[index], lat_indexes[index]] += 1
+ max_depth = numpy.max(density_map)
+
+ # create the space grids for this variable
+ space_grid = numpy.ones((grid_lon_size, grid_lat_size, max_depth), dtype=numpy.float32) * numpy.nan #TODO, dtype
+ temp_depth = numpy.zeros(space_grid_shape)
+
+ # put the variable data into the space grid
+ # FUTURE, I do not like this looping solution, figure out how to do this in native numpy ops
+ for index in range(len(data)) :
+ if not numpy.isnan(data[index]) :
+ space_grid[lon_indexes[index], lat_indexes[index], temp_depth[lon_indexes[index], lat_indexes[index]]] = data[index]
+ temp_depth[lon_indexes[index], lat_indexes[index]] += 1
+
+ return space_grid, density_map, max_depth
+
+def pack_space_grid (data_array, density_array) :
+ """
+ given a sparse array of space gridded data and a density array
+ where every slice represents the density of one "layer" in the
+ sparse array, create a well packed version of the sparse array
+ and return it
+ """
+
+ # figure out the maximum depth of the well packed data based on the density map
+ max_depth = numpy.max(numpy.sum(density_array, axis=0))
+
+ # create the final data array at the right depth
+ final_data = numpy.ones((max_depth, data_array.shape[1], data_array.shape[2]), dtype=data_array.dtype)
+ final_data = final_data * numpy.nan
+
+ LOG.debug(" original data shape: " + str(data_array.shape))
+ LOG.debug(" final data shape: " + str(final_data.shape))
+
+ # sort the sparse data into the final array
+ # FUTURE, find a way to do this without loops
+ for rows in range(data_array.shape[1]) :
+ for cols in range(data_array.shape[2]) :
+
+ # pull the data for this element
+ temp_data = data_array[:, rows, cols]
+ temp = temp_data[numpy.isfinite(temp_data)]
+
+ # put the data for this element into the final array
+ final_data[0:temp.size, rows, cols] = temp
+
+ return final_data
+
+
View
286 stg/space_time_gridding.py
@@ -0,0 +1,286 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+The main module for the Space Time Gridding application. This module handles
+the command line input from the user and coordinates the other modules.
+
+:author: Eva Schiffer (evas)
+:contact: eva.schiffer@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+
+Copyright (C) 2014 Space Science and Engineering Center (SSEC),
+ University of Wisconsin-Madison.
+"""
+__docformat__ = "restructuredtext en"
+
+import logging
+import pkg_resources
+import os
+import math
+import glob
+
+import numpy
+
+from collections import defaultdict
+
+import keoni.fbf.workspace as Workspace
+import keoni.fbf as fbf
+
+from stg.constants import *
+import stg.general_guidebook as general_guidebook
+import stg.io_manager as io_manager
+import stg.space_gridding as space_gridding
+
+# TODO, in the long run handle the dtype more flexibly
+TEMP_DATA_TYPE = numpy.dtype(numpy.float32)
+
+LOG = logging.getLogger(__name__)
+
+def get_version_string() :
+ version_num = pkg_resources.require('spacetimegrid')[0].version
+
+ return "Space Time Gridding, version " + str(version_num)
+
+def remove_file_patterns(path, *args):
+ """Remove files that were created from a previous run,
+ including temporary files, that may conflict with
+ future processing.
+
+ """
+ for pat_list in args:
+ for pat in pat_list:
+ full_pat = os.path.join(path, pat)
+ for f in glob.glob(full_pat):
+ _safe_remove(f)
+
+def _safe_remove(fn):
+ """Remove the file `fn` if you can, if not log an error message,
+ but continue on.
+
+ :Parameters:
+ fn : str
+ Filename of the file to be removed.
+ """
+ try:
+ LOG.info("Removing %s" % fn)
+ os.remove(fn)
+ except StandardError:
+ LOG.error("Could not remove %s" % fn)
+
+def main():
+ import optparse
+ usage = """
+%prog [options]
+run "%prog help" to list commands
+examples:
+
+python -m space_time_gridding
+
+"""
+
+ # set the options available to the user on the command line
+ parser = optparse.OptionParser(usage)
+
+ # message related options
+ parser.add_option('-q', '--quiet', dest="quiet",
+ action="store_true", default=False, help="only error output")
+ parser.add_option('-v', '--verbose', dest="verbose",
+ action="store_true", default=False, help="enable more informational output")
+ parser.add_option('-w', '--debug', dest="debug",
+ action="store_true", default=False, help="enable debug output")
+ parser.add_option('-n', '--version', dest='version',
+ action="store_true", default=False, help="view the STG version")
+
+ # output generation related options
+ parser.add_option('-i', '--input', dest="inputPath", type='string', default='./',
+ help="set path for the input directory")
+ parser.add_option('-o', '--output', dest="outputPath", type='string', default='./out/',
+ help="set path for the output directory")
+
+ # options related to space or time gridding
+ parser.add_option('-g', '--grid_degrees', dest="gridDegrees", type='float', default=1.0,
+ help="set the size of the output grid's cells in degrees")
+ parser.add_option('-a', '--min_scan_angle', dest="minScanAngle", type='float', default=60.0,
+ help="the minimum scan angle that will be considered useful")
+
+ # parse the uers options from the command line
+ options, args = parser.parse_args()
+
+ # set up the logging level based on the options the user selected on the command line
+ lvl = logging.WARNING
+ if options.debug: lvl = logging.DEBUG
+ elif options.verbose: lvl = logging.INFO
+ elif options.quiet: lvl = logging.ERROR
+ logging.basicConfig(level = lvl)
+
+ # display the version
+ if options.version :
+ print (get_version_string() + '\n')
+
+ commands = {}
+ prior = None
+ prior = dict(locals())
+
+ """
+ The following functions represent available menu selections
+ """
+
+ def space(*args):
+ """grid the input files in space
+ given an input directory that contains appropriate files MODIS,
+ grid them in space and put the resulting gridded date in the
+ output directory.
+
+ Note: the output directory will also be used for intermediary working
+ files.
+ """
+
+ # set up some of our input from the caller for easy access
+ desired_variables = list(args) if len(args) > 0 else [ ]
+ input_path = options.inputPath
+ output_path = options.outputPath
+ min_scan_angle = options.minScanAngle
+ grid_degrees = float(options.gridDegrees)
+
+ # determine the grid size in number of elements
+ grid_lat_size = int(math.ceil(180.0 / grid_degrees))
+ grid_lon_size = int(math.ceil(360.0 / grid_degrees))
+ space_grid_shape = (grid_lon_size, grid_lat_size) # TODO, is this the correct order?
+
+ # look through our files and figure out what variables we expect from them
+ possible_files = os.listdir(input_path)
+ expected_vars = { }
+ all_vars = set()
+ for file_name in possible_files :
+ expected_vars[file_name] = general_guidebook.get_variable_names (file_name, user_requested_names=desired_variables)
+ all_vars.update(expected_vars[file_name])
+
+ # check to make sure our intermediate file names don't exist already
+ for var_name in all_vars :
+
+ for suffix in io_manager.ALL_EXPECTED_SUFFIXES :
+ temp_name = fbf.filename(var_name + suffix, TEMP_DATA_TYPE, shape=(space_grid_shape))
+ if os.path.exists(os.path.join(output_path, temp_name)) :
+ LOG.warn ("Cannot process files because intermediate files exist in the output directory.")
+ return
+
+ # loop to deal with data from each of the files
+ for each_file in sorted(possible_files) :
+
+ full_file_path = os.path.join(input_path, each_file)
+
+ LOG.debug("Processing file: " + full_file_path)
+
+ # load the aux data
+ file_object, temp_aux_data = io_manager.load_aux_data(full_file_path,
+ min_scan_angle)
+
+ # calculate the indecies for the space grid based on the aux data
+ # (we can do this now since the lon/lat is the same for each variable in the file)
+ day_lon_index, day_lat_index, night_lon_index, night_lat_index = space_gridding.calculate_index_from_nav_data(temp_aux_data,
+ grid_degrees)
+
+ # loop to load each variable in the file and process it
+ for variable_name in expected_vars[each_file] :
+
+ LOG.debug("Processing variable: " + variable_name)
+
+ # load the variable
+ file_object, var_data = io_manager.load_variable_from_file (variable_name,
+ file_path=full_file_path,
+ file_object=file_object)
+
+ # split the variable by day/night
+ day_var_data = var_data[temp_aux_data[DAY_MASK_KEY]]
+ night_var_data = var_data[temp_aux_data[NIGHT_MASK_KEY]]
+
+ # space grid the data using the indexes we calculated earlier
+ day_space_grid, day_density_map, day_max_depth = space_gridding.space_grid_data (grid_lon_size, grid_lat_size,
+ day_var_data,
+ day_lon_index, day_lat_index)
+ night_space_grid, night_density_map, night_max_depth = space_gridding.space_grid_data (grid_lon_size, grid_lat_size,
+ night_var_data,
+ night_lon_index, night_lat_index)
+
+ # save the space grids and density info for this variable and it's density map to files
+
+ io_manager.save_data_to_file(variable_name + io_manager.DAY_TEMP_SUFFIX, space_grid_shape, output_path,
+ day_space_grid, TEMP_DATA_TYPE )
+ io_manager.save_data_to_file(variable_name + io_manager.DAY_DENSITY_TEMP_SUFFIX, space_grid_shape, output_path,
+ day_density_map, TEMP_DATA_TYPE)
+
+ io_manager.save_data_to_file(variable_name + io_manager.NIGHT_TEMP_SUFFIX, space_grid_shape, output_path,
+ night_space_grid, TEMP_DATA_TYPE )
+ io_manager.save_data_to_file(variable_name + io_manager.NIGHT_DENSITY_TEMP_SUFFIX, space_grid_shape, output_path,
+ night_density_map, TEMP_DATA_TYPE)
+
+ # make sure each file is closed when we're done with it
+ io_manager.close_file(full_file_path, file_object)
+
+ # collapse the per variable space grids to remove excess NaNs
+ for variable_name in all_vars :
+
+ LOG.debug("Packing space data for variable: " + variable_name)
+
+ # load the variable's density maps
+ var_workspace = Workspace.Workspace(dir=output_path)
+ day_var_density = var_workspace[variable_name + io_manager.DAY_DENSITY_TEMP_SUFFIX][:]
+ night_var_density = var_workspace[variable_name + io_manager.NIGHT_DENSITY_TEMP_SUFFIX][:]
+
+ # only do the day data if we have some
+ if numpy.sum(day_var_density) > 0 :
+
+ # load the sparse space grid
+ day_var_data = var_workspace[variable_name + io_manager.DAY_TEMP_SUFFIX][:]
+
+ # collapse the space grid
+ final_day_data = space_gridding.pack_space_grid(day_var_data, day_var_density)
+
+ # save the final array to an appropriately named file
+ io_manager.save_data_to_file(variable_name + io_manager.DAY_SUFFIX, space_grid_shape, output_path,
+ final_day_data, TEMP_DATA_TYPE, file_permissions="w")
+
+ else :
+ LOG.warn("No day data was found for variable " + variable_name + ". Day file will not be written.")
+
+ # only do night data if we have some
+ if numpy.sum(night_var_density) > 0 :
+
+ # load the sparse space grid
+ night_var_data = var_workspace[variable_name + io_manager.NIGHT_TEMP_SUFFIX][:]
+
+ # collapse the space grid
+ final_night_data = space_gridding.pack_space_grid(night_var_data, night_var_density)
+
+ # save the final array to an appropriately named file
+ io_manager.save_data_to_file(variable_name + io_manager.NIGHT_SUFFIX, space_grid_shape, output_path,
+ final_night_data, TEMP_DATA_TYPE, file_permissions="w")
+
+ else :
+ LOG.warn("No night data was found for variable " + variable_name + ". Night file will not be written.")
+
+ # remove the extra temporary files in the output directory
+ remove_suffixes = ["*" + p + "*" for p in io_manager.EXPECTED_TEMP_SUFFIXES]
+ remove_file_patterns(output_path, remove_suffixes)
+
+ # all the local public functions are considered part of the application, collect them up
+ commands.update(dict(x for x in locals().items() if x[0] not in prior))
+
+ # if what the user asked for is not one of our existing functions, print the help
+ if (not args) or (args[0] not in commands):
+ parser.print_help()
+ help()
+ return 9
+ else:
+ # call the function the user named, given the arguments from the command line
+ rc = locals()[args[0]](*args[1:])
+ return 0 if rc is None else rc
+
+ return 0 # it shouldn't be possible to get here any longer
+
+if __name__=='__main__':
+ sys.exit(main())
+
View
19 stg/time_gridding.py
@@ -0,0 +1,19 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+This module handles statistics across multiple times of pre-space-gridded data.
+
+:author: Eva Schiffer (evas)
+:contact: eva.schiffer@ssec.wisc.edu
+:organization: Space Science and Engineering Center (SSEC)
+:copyright: Copyright (c) 2014 University of Wisconsin SSEC. All rights reserved.
+:date: Jan 2014
+:license: GNU GPLv3
+
+Copyright (C) 2014 Space Science and Engineering Center (SSEC),
+ University of Wisconsin-Madison.
+"""
+__docformat__ = "restructuredtext en"
+
+
+

0 comments on commit 7de41a9

Please sign in to comment.
Something went wrong with that request. Please try again.