-
Michael Hiley authoredMichael Hiley authored
modis.py 12.43 KiB
#!/usr/bin/env python
# encoding: utf-8
"""
Author: Nick Bearson, nickb@ssec.wisc.edu
Copyright (c) 2012 University of Wisconsin SSEC. All rights reserved.
"""
from sensor import Sensor
from descriptors import reify
import numpy as np
from pyhdf import SD
from datetime import datetime, timedelta
import logging
log = logging.getLogger(__name__)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# FIXME: FUNCTIONS THAT MAY OR MAY NOT BELONG HERE, SHOULD THEY MOVE TO SENSOR OR A UTIL MODULE?
def _search_list_regex(alist, regex):
"""
Utility function to search a list of strings and return a list of any that match the given regular expression.
"""
import re # FIXME: move this function, and move this import when you do
retlist = []
for i in alist:
if re.match(regex, i):
retlist.append(i)
return retlist
# Think we can probably do without this...
#def _img_type(img_file):
# import re # FIXME: move this function, and move this import when you do
# if re.match('.*/?M[O|Y]D021KM', img_file):
# return '1KM'
# elif re.match('.*/?M[O|Y]D02HKM', img_file):
# return 'HKM'
# elif re.match('.*/?M[O|Y]D02QKM', img_file):
# return 'QKM'
# else:
# return None
airs_centers = np.array([11030, 12020, 13335, 13635, 13935, 14235], dtype=np.float32) / 1000 # FIXME: these are taken from modis
airs_bands = [31, 32, 33, 34, 35, 36] # FIXME: tag these as airs bands somehow
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Dictionary of band number : nanometers
b2w = {
'1' : 645,
'2' : 858.5,
'3' : 469,
'4' : 555,
'5' : 1240,
'6' : 1640,
'7' : 2130,
'8' : 412.5,
'9' : 443,
'10' : 488,
'11' : 531,
'12' : 551,
'13lo' : 667,
'13hi' : 667,
'14lo' : 678,
'14hi' : 678,
'15' : 748,
'16' : 869.5,
'17' : 905,
'18' : 936,
'19' : 940,
'20' : 3750,
'21' : 3959,
'22' : 3959,
'23' : 4050,
'24' : 4465.5,
'25' : 4515.5,
'26' : 1375,
'27' : 6715,
'28' : 7325,
'29' : 8550,
'30' : 9730,
'31' : 11030,
'32' : 12020,
'33' : 13335,
'34' : 13635,
'35' : 13935,
'36' : 14235,
}
# Add error handling here:
def nband2wave(b):
return (b2w[b] / float(1000)) # convert to micrometers (microns)
class MODIS(Sensor):
"""
MODIS sensor object to give the UWIFF implementation everything and anything it needs from the M[O|Y]D02 and M[O|Y]D03 files.
"""
SENSOR_EPOCH = datetime(1993, 1, 1)
GROUPING_REGEX = '.*M[O|Y]D03\.A([0-9]{7}\.[0-9]{4})\..*\.hdf'
GROUPING_GLOB = '.*A%s.*.hdf'
nDETECTORS = 10 # FIXME: done for 1km - does this change for hkm and qkm?
def __init__(self, tag, platform, *args, **kwargs):
# At the moment we have no extra args, but process them like this if we do...
try:
pass
# FIXME: CURRENTLY A NULL OP
# self._w = kwargs.pop('w')
except KeyError:
pass
super(MODIS, self).__init__(*args, **kwargs)
log.debug(self._file_list)
if platform == "aqua":
GEO_REGEX = '.*/?MYD03'
IMG_REGEX = '.*/?MYD02[1KM|HKM|QKM]'
CMA_REGEX = '.*/?MYD35_L2'
elif platform == "terra":
GEO_REGEX = '.*/?MOD03'
IMG_REGEX = '.*/?MOD02[1KM|HKM|QKM]'
CMA_REGEX = '.*/?MOD35_L2'
else:
log.error("MODIS platform %s not recognized! Object will not be initialized properly.")
raise RuntimeError
AIRS_REGEX = '.*/?AIRS'
self.SENSOR = "modis" # FIXME: get rid of these capitals what was i thinking
self.PLATFORM = platform
self.tag = tag
self._geo_file = _search_list_regex(self._file_list, GEO_REGEX)
# print self._geo_file
if len(self._geo_file) != 1:
raise RuntimeError, "Invalid geo list: %s" % (self._geo_file)
else:
self._geo_file = self._geo_file[0]
# FIXME: should probably separate the cloudmask more than this...
if self.tag != 'cma':
self._img_file = _search_list_regex(self._file_list, IMG_REGEX)
# print self._img_file
if len(self._img_file) != 1:
raise RuntimeError, "Invalid image list: %s" % (self._img_file)
else:
self._img_file = self._img_file[0]
if self.tag == 'cma':
self._cma_file = _search_list_regex(self._file_list, CMA_REGEX)
if len(self._cma_file) != 1:
raise RuntimeError, "Invalid cloudmask list: %s" % (self._cma_file)
else:
self._cma_file = self._cma_file[0]
self._airs_files = sorted(_search_list_regex(self._file_list, AIRS_REGEX)) # FIXME: sorting to make them line up better - should probably be linked in some other way
# self._img_type = _img_type(self._img_file)
self._refs = None
self._ref_wav = None
self._ref_nums = None
self._emis = None
self._emis_wav = None
self._emis_nums = None
@staticmethod
def date_bounds_from_fname(fname):
if "AIRS" in fname:
# FIXME: SHOULD OPEN FILE AND GET TIMES FROM INSIDE
# metastr = sd.coremetadata
# search for PRODUCTIONDATETIME
# SHORTCUT: we know we have 5 minute granules, so just make all AIRS granules overlap
start = datetime(1990, 1, 1)
end = datetime(3000, 1, 1)
return start, end
else:
DATE_REGEX = ".*A([0-9]{7}).([0-9]{4}).*.hdf"
import re
m = re.match(DATE_REGEX, fname)
if m:
d, t = m.group(1), m.group(2)
frmt = '%Y%j%H%M'
start = datetime.strptime(d+t, frmt)
end = start + timedelta(minutes=4, seconds=59)
return start, end
else:
raise RuntimeError, "can't get valid date bounds from file:\n%s\ndoesn't fit format:\n%s" % (fname, DATE_REGEX)
def _generic_info(self, ds_names, type):
ref_data = np.array([], dtype=np.float32)
ref_wavs = np.array([], dtype=np.float32)
ref_nums = []
sd = SD.SD(self._img_file)
for ds_name in ds_names:
ds = sd.select(ds_name)
valid_range = ds.valid_range
scales = np.array(getattr(ds, type + "_scales")).astype(np.float)
scales = scales.reshape(scales.shape[0], 1, 1)
offsets = np.array(getattr(ds, type + "_offsets")).astype(np.float)
offsets = offsets.reshape(offsets.shape[0], 1, 1)
vals = ds[:]
invalid_indices = (vals < valid_range[0]) | (vals > valid_range[1])
vals = (vals - offsets) * scales
vals[invalid_indices] = ds[:][invalid_indices]
bandnums = ds.band_names.split(',')
centerwaves = np.array(map(nband2wave, bandnums))
def stack(arr1, arr2):
if arr1.size: # check that arr1 is not empty
return np.concatenate((arr1, arr2))
else:
return arr2
ref_data = stack(ref_data, vals.astype(np.float32))
ref_wavs = stack(ref_wavs, centerwaves.astype(np.float32))
ref_nums.extend(bandnums)
ref_nums = np.array(ref_nums)
return (ref_data, ref_wavs, ref_nums)
@reify
def _ref_info(self):
ds_names = []
if self.tag == "1km":
ds_names.append("EV_250_Aggr1km_RefSB") # 1,2
ds_names.append("EV_500_Aggr1km_RefSB") # 3,4,5,6,7
ds_names.append("EV_1KM_RefSB") # 8,9,10,11,12,13lo,13hi,14lo,14hi,15,16,17,18,19,26
# elif self._img_type == "HKM":
# elif self._img_type == "QKM":
else:
log.error("MODIS img type %s is not currently supported" % (self.tag))
raise NotImplementedError
return self._generic_info(ds_names, "reflectance")
@reify
def _emis_info(self):
ds_names = []
if self.tag == "1km":
ds_names.append("EV_1KM_Emissive") # 20,21,22,23,24,25,27,28,29,30,31,32,33,34,35,36
# elif self._img_type == "HKM":
# elif self._img_type == "QKM":
else:
log.error("MODIS img type %s is not currently supported" % (self.tag))
raise NotImplementedError
return self._generic_info(ds_names, "radiance")
@reify
def _bt_info(self):
import modis_bright.modis_bright as mb # FIXME: packaging on this is really weird atm; mb.mb.mb to get to the actual function
emis_data, wavs, nums = self._emis_info
bt_data = np.ones(emis_data.shape) * -65535 # FIXME: should put MODIS NaN value here
for i in range(emis_data.shape[0]):
bt_data[i] = mb.modis_bright(self.PLATFORM, int(nums[i]), emis_data[i])
return (bt_data, wavs, nums)
@reify
def ReflectiveSolarBands(self):
return self._ref_info[0]
@reify
def ReflectiveBandCenters(self):
return self._ref_info[1]
@reify
def ReflectiveBandNums(self):
return self._ref_info[2]
@reify
def AirsBands(self):
import a2m
log.debug('making airs bands!')
airsshape = (len(airs_bands),) + self.Latitude.shape
allbts = np.empty(airsshape)
allbts.fill(-999.0)
allrads = np.empty(airsshape)
allrads.fill(-999.0)
for i in range(len(self._airs_files)):
airs = self._airs_files[i]
modis_geo = self._geo_file
log.debug("AirsBands processing...")
log.debug("airs file: %s" % (airs))
log.debug("modis geo file: %s" % (modis_geo))
mapping = a2m.airs2modis_transform(airs, modis_geo)
radlist = []
btlist = []
for band in airs_bands:
log.debug("adding airs band %d" % (band))
brads, bbts = a2m.airs2modis_translate(airs, band)
vrads = a2m.apply_mapping(mapping, brads)
vbts = a2m.apply_mapping(mapping, bbts)
radlist.append(vrads)
btlist.append(vbts)
rads = np.array(radlist, dtype=radlist[0].dtype)
bts = np.array(btlist, dtype=btlist[0].dtype)
# mix the new array into the total array - FIXME: REMEMBER WE'RE USING NEAREST - ADD WEIGHTING WHEN ITS AVAILABLE
def combinate(total, new):
if total == None:
log.debug("combining airs granules: first in set")
total = new
else:
assert( total.shape == new.shape ) # new should always be the same shape as total... right? catch if it's not because that could break other things
total[total == -999.0] = new[total == -999.0] # FIXME: -999.0 is just our fill - replace this with checking for larger weights?
return total
allrads = combinate(allrads, rads)
allbts = combinate(allbts, bts)
log.debug("returning from airsBands")
return allrads, allbts
def SounderEmissive(self):
return self.AirsBands[0]
def SounderBTs(self):
return self.AirsBands[1]
@reify
def SounderCenters(self):
return airs_centers
@reify
def SounderNums(self):
clist = []
for band in airs_bands:
clist.append("AIRS-%d" % (band))
return clist
@reify
def EmissiveBands(self):
return self._emis_info[0]
@reify
def EmissiveBandCenters(self):
return self._emis_info[1]
@reify
def EmissiveBandNums(self):
return self._emis_info[2]
@reify
def BrightnessTemperatureBands(self):
return self._bt_info[0]
@reify
def BrightnessTemperatureBandCenters(self):
return self._bt_info[1]
@reify
def BrightnessTemperatureBandNums(self):
return self._bt_info[2]
@reify
def Height(self):
sd = SD.SD(self._geo_file)
height = sd.select('Height')
return height[:].astype(np.float32)
@reify
def Range(self):
sd = SD.SD(self._geo_file)
r = sd.select('Range')
r = r[:] * r.scale_factor
return r.astype(np.float32)
@reify
def SensorZenith(self):
sd = SD.SD(self._geo_file)
r = sd.select('SensorZenith')
r = r[:] * r.scale_factor
return r.astype(np.float32)
@reify
def SensorAzimuth(self):
sd = SD.SD(self._geo_file)
r = sd.select('SensorAzimuth')
r = r[:] * r.scale_factor
return r.astype(np.float32)
@reify
def SolarZenith(self):
sd = SD.SD(self._geo_file)
r = sd.select('SolarZenith')
r = r[:] * r.scale_factor
return r.astype(np.float32)
@reify
def SolarAzimuth(self):
sd = SD.SD(self._geo_file)
r = sd.select('SolarAzimuth')
r = r[:] * r.scale_factor
return r.astype(np.float32)
@reify
def Latitude(self):
sd = SD.SD(self._geo_file)
r = sd.select('Latitude')
return r[:]
@reify
def Longitude(self):
sd = SD.SD(self._geo_file)
r = sd.select('Longitude')
return r[:]
@reify
def ScanStartTime(self):
sd = SD.SD(self._geo_file)
r = sd.select('SD start time')
return r[:]
@reify
def LandSeaMask_compare(self):
print "GETTING LANDSEAMASK FROM FILE"
sd = SD.SD(self._geo_file)
r = sd.select('Land/SeaMask')
return r[:]
def Cloudmask(self, qf):
sd = SD.SD(self._cma_file)
r = sd.select('Cloud_Mask')
i = qf - 1 # qf is 1-6, we want 0-based
return r[i]
@reify
def CMQA(self):
sd = SD.SD(self._cma_file)
r = sd.select('Quality_Assurance')
return r[:]