diff --git a/COTconvert.py b/COTconvert.py
deleted file mode 100644
index 1d073c9f8d3249604ef25c545d5a2480e55ded28..0000000000000000000000000000000000000000
--- a/COTconvert.py
+++ /dev/null
@@ -1,103 +0,0 @@
-#!/usr/bin/env python
-
-### System libraries
-import time, string, sys
-from os import path
-
-import numpy as np
-from bisect import bisect,bisect_right,bisect_left
-
-def binary_search_SP(cpsval,size,cpsArr):
- status = 0
- if ((cpsval > cpsArr[0]) and (cpsval < cpsArr[-1])):
- #-- This is for unit offset
- #l_index = 0
- #u_index = size+1
- #-- This is for zero offset
- l_index = -1
- u_index = size
- while ((u_index - l_index) > 1):
- search_index = (u_index + l_index)/2
- if (cpsval <= cpsArr[search_index]):
- u_index = search_index
- else :
- l_index = search_index
-
- elif (cpsval <= cpsArr[0]):
- print "cpsval <= minimum in LUT"
- l_index = 0
- u_index = 1
- status = 1
- elif (cpsval >= cpsArr[-1]):
- print "cpsval >= maximum in LUT"
- l_index = size-2
- u_index = size-1
- status = 1
- else :
- print "Some situation not captured, passing!!!"
- pass
-
- return [u_index, l_index, status]
-
-def COTconvert(nCOT,COTlut,CPS,COTvis):
-
- cpsArr = COTlut[:,0]
- visIrEfficRatio = COTlut[:,1]
-
- # If CPS exceeds largest value in LUT, assign to COTir
- # the last COT value in LUT
-
- if ( CPS > cpsArr[-1]) :
- s = visIrEfficRatio[-1]
- else :
- u_index,l_index,status = binary_search_SP(CPS, nCOT, cpsArr)
-
- print "Finished binary search..."
- print "l_index = ",l_index
- print "u_index = ",u_index
-
- m = (visIrEfficRatio[u_index]-visIrEfficRatio[l_index])/(cpsArr[u_index]-cpsArr[l_index])
- b = visIrEfficRatio[u_index]-m*cpsArr[u_index]
- s = m*CPS+b
-
- print "cpsArr[l_index] = ",cpsArr[l_index]
- print "cpsArr[u_index] = ",cpsArr[u_index]
- print "visIrEfficRatio[l_index] = ",visIrEfficRatio[l_index]
- print "visIrEfficRatio[u_index] = ",visIrEfficRatio[u_index]
- print "s = ",s
-
- # Compute the IR COT
- COTir = COTvis/s
- return COTir
-
-def readCOTfile(file):
- results = []
- file=str(file)
- f = open(file,'r')
-
- lines = f.readlines()[1:]
- f.close()
- for l in lines:
- fields = l.split()
- CPS = float(fields[0])
- COT = float(fields[1])
- all = [CPS,COT]
- results.append(all)
-
- results = np.array(results)
- nCOT = np.shape(results)[0]
- return results,nCOT
-
-def main(cot,cps):
- file = '/data/geoffc/MODIS_NPOESS/trunk/LEOCAT/data/VIIRS/LUT/CTP/COTlut.db'
- COTlut,nCOT, = readCOTfile(file)
- COTir = COTconvert(nCOT,COTlut,cps,cot)
- print "COTir = ",COTir
-
-if __name__=='__main__':
- cot = float(sys.argv[1])
- cps = float(sys.argv[2])
- print "Input cloud optical thickness: ",cot
- print "Input cloud particle size: ",cps,"microns"
-
- main(cot,cps)
diff --git a/LandWaterMask.py b/LandWaterMask.py
new file mode 100644
index 0000000000000000000000000000000000000000..f7cff86f97de15dcf53992ab19ca96b6a97523ef
--- /dev/null
+++ b/LandWaterMask.py
@@ -0,0 +1,527 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+LandWaterMask.py
+
+ * DESCRIPTION: Class to granulate the DEM Land Water Mask data product
+
+Created by Geoff Cureton on 2013-03-05.
+Copyright (c) 2013 University of Wisconsin SSEC. All rights reserved.
+"""
+
+file_Date = '$Date: 2014-09-08 17:01:37 -0500 (Mon, 08 Sep 2014) $'
+file_Revision = '$Revision: 2224 $'
+file_Author = '$Author: geoffc $'
+file_HeadURL = '$HeadURL: https://svn.ssec.wisc.edu/repos/jpss_adl/trunk/scripts/edr/GridIP/LandWaterMask.py $'
+file_Id = '$Id: LandWaterMask.py 2224 2014-09-08 22:01:37Z geoffc $'
+
+__author__ = 'G.P. Cureton <geoff.cureton@ssec.wisc.edu>'
+__version__ = '$Id: LandWaterMask.py 2224 2014-09-08 22:01:37Z geoffc $'
+__docformat__ = 'Epytext'
+
+
+
+import os, sys, logging, traceback
+from os import path,uname,environ
+import string
+import re
+import uuid
+from glob import glob
+from time import time
+from datetime import datetime,timedelta
+
+from scipy import round_
+
+import numpy as np
+from numpy import ma
+import copy
+from bisect import bisect_left,bisect_right
+
+import ctypes
+from numpy.ctypeslib import ndpointer
+
+import tables as pytables
+from tables import exceptions as pyEx
+
+import ViirsData
+
+# skim and convert routines for reading .asc metadata fields of interest
+import adl_blob2 as adl_blob
+import adl_asc
+from adl_asc import skim_dir, contiguous_granule_groups, granule_groups_contain, effective_anc_contains,eliminate_duplicates,_is_contiguous, RDR_REQUIRED_KEYS, POLARWANDER_REQUIRED_KEYS
+from adl_common import ADL_HOME, CSPP_RT_HOME, CSPP_RT_ANC_HOME
+
+# every module should have a LOG object
+try :
+ sourcename= file_Id.split(" ")
+ LOG = logging.getLogger(sourcename[1])
+except :
+ LOG = logging.getLogger('LandWaterMask')
+
+from Utils import getURID, getAscLine, getAscStructs, findDatelineCrossings, shipOutToFile
+from Utils import index, find_lt, find_le, find_gt, find_ge
+from Utils import plotArr
+
+
+class LandWaterMask() :
+
+ def __init__(self,inDir=None, sdrEndian=None, ancEndian=None):
+ self.collectionShortName = 'VIIRS-GridIP-VIIRS-Lwm-Mod-Gran'
+ self.xmlName = 'VIIRS_GRIDIP_VIIRS_LWM_MOD_GRAN.xml'
+ self.blobDatasetName = 'landWaterMask'
+ self.dataType = 'uint8'
+ self.sourceType = 'DEM'
+ self.sourceList = ['']
+ self.trimObj = ViirsData.ViirsTrimTable()
+
+ if inDir is None :
+ self.inDir = path.abspath(path.curdir)
+ else :
+ self.inDir = inDir
+
+ if sdrEndian is None :
+ self.sdrEndian = adl_blob.LITTLE_ENDIAN
+ else :
+ self.sdrEndian = sdrEndian
+
+ if ancEndian is None :
+ self.ancEndian = adl_blob.LITTLE_ENDIAN
+ else :
+ self.ancEndian = ancEndian
+
+ # Digital Elevation Model (DEM) land sea mask types
+ self.DEM_list = ['DEM_SHALLOW_OCEAN','DEM_LAND','DEM_COASTLINE',
+ 'DEM_SHALLOW_INLAND_WATER','DEM_EPHEMERAL_WATER',
+ 'DEM_DEEP_INLAND_WATER','DEM_MOD_CONT_OCEAN','DEM_DEEP_OCEAN']
+ self.DEM_dict = {
+ 'DEM_SHALLOW_OCEAN' : 0,
+ 'DEM_LAND' : 1,
+ 'DEM_COASTLINE' : 2,
+ 'DEM_SHALLOW_INLAND_WATER' : 3,
+ 'DEM_EPHEMERAL_WATER' : 4,
+ 'DEM_DEEP_INLAND_WATER' : 5,
+ 'DEM_MOD_CONT_OCEAN' : 6,
+ 'DEM_DEEP_OCEAN' : 7
+ }
+
+
+ def setGeolocationInfo(self,dicts):
+ '''
+ Populate this class instance with the geolocation data for a single granule
+ '''
+ # Set some environment variables and paths
+ ANC_SCRIPTS_PATH = path.join(CSPP_RT_HOME,'viirs')
+ ADL_ASC_TEMPLATES = path.join(ANC_SCRIPTS_PATH,'asc_templates')
+
+ # Collect some data from the geolocation dictionary
+ self.geoDict = dicts
+ URID = dicts['URID']
+ geo_Collection_ShortName = dicts['N_Collection_Short_Name']
+ N_Granule_ID = dicts['N_Granule_ID']
+ ObservedStartTimeObj = dicts['ObservedStartTime']
+ geoFiles = glob('%s/%s*' % (self.inDir,URID))
+ geoFiles.sort()
+
+ LOG.debug("###########################")
+ LOG.debug(" Geolocation Information ")
+ LOG.debug("###########################")
+ LOG.debug("N_Granule_ID : %r" % (N_Granule_ID))
+ LOG.debug("ObservedStartTime : %s" % (ObservedStartTimeObj.__str__()))
+ LOG.debug("N_Collection_Short_Name : %s" %(geo_Collection_ShortName))
+ LOG.debug("URID : %r" % (URID))
+ LOG.debug("geoFiles : %r" % (geoFiles))
+ LOG.debug("###########################")
+
+ # Do we have terrain corrected geolocation?
+ terrainCorrectedGeo = True if 'GEO-TC' in geo_Collection_ShortName else False
+
+ # Do we have long or short style geolocation field names?
+ if (geo_Collection_ShortName=='VIIRS-MOD-GEO-TC' or geo_Collection_ShortName=='VIIRS-MOD-RGEO') :
+ longFormGeoNames = True
+ LOG.debug("We have long form geolocation names")
+ elif (geo_Collection_ShortName=='VIIRS-MOD-GEO' or geo_Collection_ShortName=='VIIRS-MOD-RGEO-TC') :
+ LOG.debug("We have short form geolocation names")
+ longFormGeoNames = False
+ else :
+ LOG.error("Invalid geolocation shortname: %s" %(geo_Collection_ShortName))
+ return -1
+
+ # Get the geolocation xml file
+
+ geoXmlFile = "%s.xml" % (string.replace(geo_Collection_ShortName,'-','_'))
+ geoXmlFile = path.join(ADL_HOME,'xml/VIIRS',geoXmlFile)
+ if path.exists(geoXmlFile):
+ LOG.debug("We are using for %s: %s,%s" %(geo_Collection_ShortName,geoXmlFile,geoFiles[0]))
+
+ # Open the geolocation blob and get the latitude and longitude
+
+ endian = self.sdrEndian
+
+ geoBlobObj = adl_blob.map(geoXmlFile,geoFiles[0], endian=endian)
+
+ # Get scan_mode to find any bad scans
+
+ scanMode = geoBlobObj.scan_mode[:]
+ badScanIdx = np.where(scanMode==254)[0]
+ LOG.debug("Bad Scans: %r" % (badScanIdx))
+
+ # Detemine the min, max and range of the latitude and longitude,
+ # taking care to exclude any fill values.
+
+ if longFormGeoNames :
+ if endian==adl_blob.BIG_ENDIAN:
+ latitude = getattr(geoBlobObj,'latitude').byteswap()
+ longitude = getattr(geoBlobObj,'longitude').byteswap()
+ latitude = latitude.astype('float')
+ longitude = longitude.astype('float')
+ else:
+ latitude = getattr(geoBlobObj,'latitude').astype('float')
+ longitude = getattr(geoBlobObj,'longitude').astype('float')
+ else :
+ latitude = getattr(geoBlobObj,'lat').astype('float')
+ longitude = getattr(geoBlobObj,'lon').astype('float')
+
+ latitude = ma.masked_less(latitude,-800.)
+ latMin,latMax = np.min(latitude),np.max(latitude)
+ latRange = latMax-latMin
+
+ longitude = ma.masked_less(longitude,-800.)
+ lonMin,lonMax = np.min(longitude),np.max(longitude)
+ lonRange = lonMax-lonMin
+
+ LOG.debug("min,max,range of latitide: %f %f %f" % (latMin,latMax,latRange))
+ LOG.debug("min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+
+ # Determine the latitude and longitude fill masks, so we can restore the
+ # fill values after we have scaled...
+
+ latMask = latitude.mask
+ lonMask = longitude.mask
+
+ # Check if the geolocation is in radians, convert to degrees
+ if 'RGEO' in geo_Collection_ShortName :
+ LOG.debug("Geolocation is in radians, convert to degrees...")
+ latitude = np.degrees(latitude)
+ longitude = np.degrees(longitude)
+
+ latMin,latMax = np.min(latitude),np.max(latitude)
+ latRange = latMax-latMin
+ lonMin,lonMax = np.min(longitude),np.max(longitude)
+ lonRange = lonMax-lonMin
+
+ LOG.debug("New min,max,range of latitude: %f %f %f" % (latMin,latMax,latRange))
+ LOG.debug("New min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+
+ # Restore fill values to masked pixels in geolocation
+
+ geoFillValue = self.trimObj.sdrTypeFill['VDNE_FLOAT64_FILL'][latitude.dtype.name]
+ latitude = ma.array(latitude,mask=latMask,fill_value=geoFillValue)
+ self.latitude = latitude.filled()
+
+ geoFillValue = self.trimObj.sdrTypeFill['VDNE_FLOAT64_FILL'][longitude.dtype.name]
+ longitude = ma.array(longitude,mask=lonMask,fill_value=geoFillValue)
+ self.longitude = longitude.filled()
+
+ # Shift the longitudes to be between -180 and 180 degrees
+ if lonMax > 180. :
+ LOG.debug("\nFinal min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+ dateLineIdx = np.where(longitude>180.)
+ LOG.debug("dateLineIdx = %r" % (dateLineIdx))
+ longitude[dateLineIdx] -= 360.
+ lonMax = np.max(ma.array(longitude,mask=lonMask))
+ lonMin = np.min(ma.array(longitude,mask=lonMask))
+ lonRange = lonMax-lonMin
+ LOG.debug("\nFinal min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+
+ # Record the corners, taking care to exclude any bad scans...
+ nDetectors = 16
+ firstGoodScan = np.where(scanMode<=2)[0][0]
+ lastGoodScan = np.where(scanMode<=2)[0][-1]
+ firstGoodRow = firstGoodScan * nDetectors
+ lastGoodRow = lastGoodScan * nDetectors + nDetectors - 1
+
+ latCrnList = [latitude[firstGoodRow,0],latitude[firstGoodRow,-1],latitude[lastGoodRow,0],latitude[lastGoodRow,-1]]
+ lonCrnList = [longitude[firstGoodRow,0],longitude[firstGoodRow,-1],longitude[lastGoodRow,0],longitude[lastGoodRow,-1]]
+
+ # Check for dateline/pole crossings
+ num180Crossings = findDatelineCrossings(latCrnList,lonCrnList)
+ LOG.debug("We have %d dateline crossings."%(num180Crossings))
+
+ # Copy the geolocation information to the class object
+ self.latMin = latMin
+ self.latMax = latMax
+ self.latRange = latRange
+ self.lonMin = lonMin
+ self.lonMax = lonMax
+ self.lonRange = lonRange
+ self.scanMode = scanMode
+ self.latitude = latitude
+ self.longitude = longitude
+ self.latCrnList = latCrnList
+ self.lonCrnList = lonCrnList
+ self.num180Crossings = num180Crossings
+
+ # Parse the geolocation asc file to get struct information which will be
+ # written to the ancillary asc files
+
+ geoAscFileName = path.join(self.inDir,URID+".asc")
+ LOG.debug("\nOpening %s..." % (geoAscFileName))
+
+ geoAscFile = open(geoAscFileName,'rt')
+
+ self.ObservedDateTimeStr = getAscLine(geoAscFile,"ObservedDateTime")
+ #self.RangeDateTimeStr = self.ObservedDateTimeStr
+ self.RangeDateTimeStr = getAscLine(geoAscFile,"ObservedDateTime")
+ self.RangeDateTimeStr = string.replace(self.RangeDateTimeStr,"ObservedDateTime","RangeDateTime")
+ self.GRingLatitudeStr = getAscStructs(geoAscFile,"GRingLatitude",12)
+ self.GRingLongitudeStr = getAscStructs(geoAscFile,"GRingLongitude",12)
+
+ self.North_Bounding_Coordinate_Str = getAscLine(geoAscFile,"North_Bounding_Coordinate")
+ self.South_Bounding_Coordinate_Str = getAscLine(geoAscFile,"South_Bounding_Coordinate")
+ self.East_Bounding_Coordinate_Str = getAscLine(geoAscFile,"East_Bounding_Coordinate")
+ self.West_Bounding_Coordinate_Str = getAscLine(geoAscFile,"West_Bounding_Coordinate")
+
+ geoAscFile.close()
+
+
+ def subset(self):
+ '''Subsets the LSM dataset to cover the required geolocation range.'''
+
+ # Get the subset of DEM global dataset.
+
+ DEM_dLat = 30.*(1./3600.)
+ DEM_dLon = 30.*(1./3600.)
+
+ DEM_fileName = path.join(CSPP_RT_ANC_HOME,'LSM/dem30ARC_Global_LandWater_uncompressed.h5')
+ self.sourceList.append(path.basename(DEM_fileName))
+
+ try :
+ # TODO : Use original HDF4 file which contains elevation and LWM.
+ DEMobj = pytables.openFile(DEM_fileName,'r')
+ DEM_node = DEMobj.getNode('/demGRID/Data Fields/LandWater')
+ except Exception, err :
+ LOG.exception("%s"%(err))
+ LOG.exception("Problem opening DEM file (%s), aborting."%(DEM_fileName))
+ sys.exit(1)
+
+ try :
+ DEM_gridLats = -1. * (np.arange(21600.) * DEM_dLat - 90.)
+ DEM_gridLons = np.arange(43200.) * DEM_dLon - 180.
+
+ LOG.debug("min,max DEM Grid Latitude values : %f,%f"%(DEM_gridLats[0],DEM_gridLats[-1]))
+ LOG.debug("min,max DEM Grid Longitude values : %f,%f"%(DEM_gridLons[0],DEM_gridLons[-1]))
+
+ latMin = self.latMin
+ latMax = self.latMax
+ lonMin = self.lonMin
+ lonMax = self.lonMax
+
+ DEM_latMask = np.equal((DEM_gridLats<(latMax+DEM_dLat)),(DEM_gridLats>(latMin-DEM_dLat)))
+ DEM_lonMask = np.equal((DEM_gridLons<(lonMax+DEM_dLon)),(DEM_gridLons>(lonMin-DEM_dLon)))
+
+ DEM_latIdx = np.where(DEM_latMask==True)[0]
+ DEM_lonIdx = np.where(DEM_lonMask==True)[0]
+
+ DEM_latMinIdx = DEM_latIdx[0]
+ DEM_latMaxIdx = DEM_latIdx[-1]
+ DEM_lonMinIdx = DEM_lonIdx[0]
+ DEM_lonMaxIdx = DEM_lonIdx[-1]
+
+ LOG.debug("DEM_latMinIdx = %d" % (DEM_latMinIdx))
+ LOG.debug("DEM_latMaxIdx = %d" % (DEM_latMaxIdx))
+ LOG.debug("DEM_lonMinIdx = %d" % (DEM_lonMinIdx))
+ LOG.debug("DEM_lonMaxIdx = %d" % (DEM_lonMaxIdx))
+
+ lat_subset = DEM_gridLats[DEM_latMinIdx:DEM_latMaxIdx+1]
+ self.gridLat = lat_subset
+
+ if self.num180Crossings == 2 :
+
+ # We have a dateline crossing, so subset the positude and negative
+ # longitude grids and sandwich them together.
+ posLonCrn = np.min(ma.masked_less_equal(np.array(self.lonCrnList),0.))
+ negLonCrn = np.max(ma.masked_outside(np.array(self.lonCrnList),-800.,0.))
+ posIdx = index(DEM_gridLons,find_lt(DEM_gridLons,posLonCrn))
+ negIdx = index(DEM_gridLons,find_gt(DEM_gridLons,negLonCrn))
+
+ posLons_subset = DEM_gridLons[posIdx:]
+ negLons_subset = DEM_gridLons[:negIdx]
+ lon_subset = np.concatenate((posLons_subset,negLons_subset))
+
+ # Do the same with the DEM data
+ posBlock = DEM_node[DEM_latMinIdx:DEM_latMaxIdx+1,posIdx:]
+ negBlock = DEM_node[DEM_latMinIdx:DEM_latMaxIdx+1,:negIdx]
+ DEM_subset = np.concatenate((posBlock,negBlock),axis=1)
+
+ else :
+
+ DEM_subset = DEM_node[DEM_latMinIdx:DEM_latMaxIdx+1,DEM_lonMinIdx:DEM_lonMaxIdx+1]
+ lon_subset = DEM_gridLons[DEM_lonMinIdx:DEM_lonMaxIdx+1]
+
+ self.gridLon = lon_subset
+
+ # Copy DEM data to the GridIP object
+ self.gridData = DEM_subset.astype(self.dataType)
+
+ DEM_node.close()
+ DEMobj.close()
+
+ except Exception, err :
+
+ LOG.debug("EXCEPTION: %s" % (err))
+
+ DEM_node.close()
+ DEMobj.close()
+
+
+ def _grid2Gran(self, dataLat, dataLon, gridData, gridLat, gridLon):
+ '''Granulates a gridded dataset using an input geolocation'''
+
+ nData = np.int64(dataLat.size)
+ gridRows = np.int32(gridLat.shape[0])
+ gridCols = np.int32(gridLat.shape[1])
+
+ data = np.ones(np.shape(dataLat),dtype=np.float64)* -999.9
+ dataIdx = np.ones(np.shape(dataLat),dtype=np.int64) * -99999
+
+ ANC_SCRIPTS_PATH = path.join(CSPP_RT_HOME,'viirs')
+
+ libFile = path.join(ANC_SCRIPTS_PATH,'libgriddingAndGranulation.so.1.0.1')
+ LOG.debug("Gridding and granulation library file: %s" % (libFile))
+ lib = ctypes.cdll.LoadLibrary(libFile)
+ grid2gran = lib.grid2gran_nearest
+ grid2gran.restype = None
+ grid2gran.argtypes = [
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ctypes.c_int64,
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_int64,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ctypes.c_int32,
+ ctypes.c_int32
+ ]
+
+ '''
+ int snapGrid_ctypes(double *lat,
+ double *lon,
+ double *data,
+ long nData,
+ double *gridLat,
+ double *gridLon,
+ double *gridData,
+ long *gridDataIdx,
+ int nGridRows,
+ int nGridCols
+ )
+ '''
+
+ LOG.debug("Calling C routine grid2gran()...")
+
+ retVal = grid2gran(dataLat,
+ dataLon,
+ data,
+ nData,
+ gridLat,
+ gridLon,
+ gridData,
+ dataIdx,
+ gridRows,
+ gridCols)
+
+ LOG.debug("Returning from C routine grid2gran()")
+
+
+ return data,dataIdx
+
+
+ def granulate(self,GridIP_objects):
+ '''
+ Granulates the GridIP DEM files.
+ '''
+
+ # Generate the lat and lon grids, and flip them and the data over latitude
+ gridLon,gridLat = np.meshgrid(self.gridLon,self.gridLat[::-1])
+ gridData = self.gridData[::-1,:]
+
+ latitude = self.latitude
+ longitude = self.longitude
+
+ # If we have a dateline crossing, remove the longitude discontinuity
+ # by adding 360 degrees to the negative longitudes.
+ if self.num180Crossings == 2 :
+ gridLonNegIdx = np.where(gridLon < 0.)
+ gridLon[gridLonNegIdx] += 360.
+ longitudeNegIdx = np.where(longitude < 0.)
+ longitude[longitudeNegIdx] += 360.
+
+ LOG.info("Granulating %s ..." % (self.collectionShortName))
+ LOG.debug("latitide,longitude shapes: %s, %s"%(str(latitude.shape) , str(longitude.shape)))
+ LOG.debug("gridData.shape = %s" % (str(gridData.shape)))
+ LOG.debug("gridLat.shape = %s" % (str(gridLat.shape)))
+ LOG.debug("gridLon.shape = %s" % (str(gridLon.shape)))
+
+ LOG.debug("min of gridData = %r"%(np.min(gridData)))
+ LOG.debug("max of gridData = %r"%(np.max(gridData)))
+
+ t1 = time()
+ data,dataIdx = self._grid2Gran(np.ravel(latitude),
+ np.ravel(longitude),
+ gridData.astype(np.float64),
+ gridLat.astype(np.float64),
+ gridLon.astype(np.float64))
+ t2 = time()
+ elapsedTime = t2-t1
+ LOG.info("Granulation took %f seconds for %d points" % (elapsedTime,latitude.size))
+
+ data = data.reshape(latitude.shape)
+ dataIdx = dataIdx.reshape(latitude.shape)
+
+ # Convert granulated data back to original type...
+ data = data.astype(self.dataType)
+
+ # Convert any "inland water" to "sea water"
+ shallowInlandWaterValue = self.DEM_dict['DEM_SHALLOW_INLAND_WATER']
+ shallowOceanValue = self.DEM_dict['DEM_SHALLOW_OCEAN']
+ deepInlandWaterValue = self.DEM_dict['DEM_DEEP_INLAND_WATER']
+ deepOceanValue = self.DEM_dict['DEM_DEEP_OCEAN']
+
+ shallowInlandWaterMask = ma.masked_equal(data,shallowInlandWaterValue).mask
+ shallowOceanMask = ma.masked_equal(data,shallowOceanValue).mask
+ deepInlandWaterMask = ma.masked_equal(data,deepInlandWaterValue).mask
+
+ totalWaterMask = shallowInlandWaterMask #+ shallowOceanMask + deepInlandWaterMask
+
+ data = ma.array(data,mask=totalWaterMask,fill_value=deepOceanValue)
+ data = data.filled()
+
+
+ LOG.debug("Shape of granulated %s data is %s" % (self.collectionShortName,np.shape(data)))
+ LOG.debug("Shape of granulated %s dataIdx is %s" % (self.collectionShortName,np.shape(dataIdx)))
+
+ # Explicitly restore geolocation fill to the granulated data...
+ fillMask = ma.masked_less(self.latitude,-800.).mask
+ fillValue = self.trimObj.sdrTypeFill['MISS_FILL'][self.dataType]
+ data = ma.array(data,mask=fillMask,fill_value=fillValue)
+ self.data = data.filled()
+
+ # Moderate resolution trim table arrays. These are
+ # bool arrays, and the trim pixels are set to True.
+ #modTrimMask = self.trimObj.createModTrimArray(nscans=48,trimType=bool)
+
+ # Fill the required pixel trim rows in the granulated GridIP data with
+ # the ONBOARD_PT_FILL value for the correct data type
+
+ #fillValue = self.trimObj.sdrTypeFill['ONBOARD_PT_FILL'][self.dataType]
+ #data = ma.array(data,mask=modTrimMask,fill_value=fillValue)
+ #self.data = data.filled()
+
+
+ def shipOutToFile(self):
+ ''' Pass the current class instance to this Utils method to generate
+ a blob/asc file pair from the input ancillary data object.'''
+
+ #shipOutToFile(self)
+ pass
diff --git a/PrecipWater.py b/PrecipWater.py
new file mode 100644
index 0000000000000000000000000000000000000000..c34c81fc6f651c128f440ef38c914c8b19df5349
--- /dev/null
+++ b/PrecipWater.py
@@ -0,0 +1,434 @@
+#!/usr/bin/env python
+# encoding: utf-8
+"""
+PrecipWater.py
+
+ * DESCRIPTION: Class to granulate the ViirsAncPrecipWater data product
+
+Created by Geoff Cureton on 2013-02-25.
+Copyright (c) 2011 University of Wisconsin SSEC. All rights reserved.
+"""
+
+file_Date = '$Date: 2014-09-19 15:13:47 -0500 (Fri, 19 Sep 2014) $'
+file_Revision = '$Revision: 2235 $'
+file_Author = '$Author: geoffc $'
+file_HeadURL = '$HeadURL: https://svn.ssec.wisc.edu/repos/jpss_adl/trunk/scripts/edr/ANC/PrecipWater.py $'
+file_Id = '$Id: PrecipWater.py 2235 2014-09-19 20:13:47Z geoffc $'
+
+__author__ = 'G.P. Cureton <geoff.cureton@ssec.wisc.edu>'
+__version__ = '$Id: PrecipWater.py 2235 2014-09-19 20:13:47Z geoffc $'
+__docformat__ = 'Epytext'
+
+
+
+import os, sys, logging, traceback
+from os import path,uname,environ
+import string
+import re
+import uuid
+from glob import glob
+from time import time
+from datetime import datetime,timedelta
+
+from scipy import round_
+
+import numpy as np
+from numpy import ma
+import copy
+from bisect import bisect_left,bisect_right
+
+import ctypes
+from numpy.ctypeslib import ndpointer
+
+import ViirsData
+
+from NCEPtoBlob import NCEPclass
+
+# skim and convert routines for reading .asc metadata fields of interest
+import adl_blob2 as adl_blob
+import adl_asc
+from adl_asc import skim_dir, contiguous_granule_groups, granule_groups_contain, \
+ effective_anc_contains,eliminate_duplicates,_is_contiguous, \
+ RDR_REQUIRED_KEYS, POLARWANDER_REQUIRED_KEYS
+from adl_common import ADL_HOME, CSPP_RT_HOME, CSPP_RT_ANC_PATH, \
+ CSPP_RT_ANC_CACHE_DIR, COMMON_LOG_CHECK_TABLE
+
+# every module should have a LOG object
+try :
+ sourcename= file_Id.split(" ")
+ LOG = logging.getLogger(sourcename[1])
+except :
+ LOG = logging.getLogger('PrecipWater')
+
+from Utils import getURID, getAscLine, getAscStructs, findDatelineCrossings, \
+ shipOutToFile
+
+
+class PrecipWater() :
+
+ def __init__(self,inDir=None, sdrEndian=None, ancEndian=None):
+ self.collectionShortName = 'VIIRS-ANC-Preci-Wtr-Mod-Gran'
+ self.xmlName = 'VIIRS_ANC_PRECI_WTR_MOD_GRAN.xml'
+ self.blobDatasetName = 'totalPrecipitableWater'
+ self.dataType = 'float32'
+ self.sourceType = 'NCEP_ANC_Int'
+ self.sourceList = ['']
+ self.trimObj = ViirsData.ViirsTrimTable()
+
+ if inDir is None :
+ self.inDir = path.abspath(path.curdir)
+ else :
+ self.inDir = inDir
+
+ if sdrEndian is None :
+ self.sdrEndian = adl_blob.LITTLE_ENDIAN
+ else :
+ self.sdrEndian = sdrEndian
+
+ if ancEndian is None :
+ self.ancEndian = adl_blob.LITTLE_ENDIAN
+ else :
+ self.ancEndian = ancEndian
+
+
+ def ingest(self,ancBlob=None):
+ '''
+ Ingest the ancillary dataset.
+ '''
+ dates = []
+ ncepBlobFiles = []
+ for gridBlobStruct in ancBlob:
+ timeObj = gridBlobStruct[0]
+ ncepBlobFile = gridBlobStruct[1]
+ LOG.debug("VIIRS-ANC-Preci-Wtr-Mod-Gran %s --> %s" % \
+ (ncepBlobFile,timeObj.strftime("%Y-%m-%d %H:%M:%S:%f")))
+ dates.append(timeObj)
+ ncepBlobFiles.append(ncepBlobFile)
+
+ self.date_0 = dates[0]
+ self.date_1 = dates[1]
+
+ LOG.debug("Minimum NCEP date is: %s" %(self.date_0.strftime("%Y-%m-%d %H:%M:%S:%f")))
+ LOG.debug("Maximum NCEP date is: %s" %(self.date_1.strftime("%Y-%m-%d %H:%M:%S:%f")))
+
+ ncepBlobFile_0 = ncepBlobFiles[0]
+ ncepBlobFile_1 = ncepBlobFiles[1]
+
+ self.gridData_0 = getattr(ncepBlobFile_0,self.blobDatasetName).astype(self.dataType)
+ self.gridData_1 = getattr(ncepBlobFile_1,self.blobDatasetName).astype(self.dataType)
+
+
+ def setGeolocationInfo(self,dicts):
+ '''
+ Populate this class instance with the geolocation data for a single granule
+ '''
+ # Set some environment variables and paths
+ ANC_SCRIPTS_PATH = path.join(CSPP_RT_HOME,'viirs')
+ ADL_ASC_TEMPLATES = path.join(ANC_SCRIPTS_PATH,'asc_templates')
+
+ # Collect some data from the geolocation dictionary
+ self.geoDict = dicts
+ URID = dicts['URID']
+ geo_Collection_ShortName = dicts['N_Collection_Short_Name']
+ N_Granule_ID = dicts['N_Granule_ID']
+ ObservedStartTimeObj = dicts['ObservedStartTime']
+ geoFiles = glob('%s/%s*' % (self.inDir,URID))
+ geoFiles.sort()
+
+ LOG.debug("###########################")
+ LOG.debug(" Geolocation Information ")
+ LOG.debug("###########################")
+ LOG.debug("N_Granule_ID : %r" % (N_Granule_ID))
+ LOG.debug("ObservedStartTime : %s" % (ObservedStartTimeObj.__str__()))
+ LOG.debug("N_Collection_Short_Name : %s" %(geo_Collection_ShortName))
+ LOG.debug("URID : %r" % (URID))
+ LOG.debug("geoFiles : %r" % (geoFiles))
+ LOG.debug("###########################")
+
+ timeDelta = (self.date_1 - self.date_0).total_seconds()
+ LOG.debug("timeDelta is %r seconds" %(timeDelta))
+
+ timePrime = (ObservedStartTimeObj - self.date_0).total_seconds()
+ LOG.debug("timePrime is %r seconds (%f percent along time interval)" % \
+ (timePrime,(timePrime/timeDelta)*100.))
+
+ delta_gridData = self.gridData_1 - self.gridData_0
+
+ self.gridData = (delta_gridData/timeDelta) * timePrime + self.gridData_0
+
+ gridData_0_avg = np.average(self.gridData_0)
+ gridData_1_avg = np.average(self.gridData_1)
+ gridData_avg = np.average(self.gridData)
+ LOG.debug("average(gridData_0) = %f" %(np.average(self.gridData_0)))
+ LOG.debug("average(gridData_1) = %f" %(np.average(self.gridData_1)))
+ LOG.debug("average(gridData) = %f" %(np.average(self.gridData)))
+
+ # Do we have terrain corrected geolocation?
+ terrainCorrectedGeo = True if 'GEO-TC' in geo_Collection_ShortName else False
+
+ # Do we have long or short style geolocation field names?
+ if (geo_Collection_ShortName=='VIIRS-MOD-GEO-TC' or geo_Collection_ShortName=='VIIRS-MOD-RGEO') :
+ longFormGeoNames = True
+ LOG.debug("We have long form geolocation names")
+ elif (geo_Collection_ShortName=='VIIRS-MOD-GEO' or geo_Collection_ShortName=='VIIRS-MOD-RGEO-TC') :
+ LOG.debug("We have short form geolocation names")
+ longFormGeoNames = False
+ else :
+ LOG.error("Invalid geolocation shortname: %s" %(geo_Collection_ShortName))
+ return -1
+
+ # Get the geolocation xml file
+
+ geoXmlFile = "%s.xml" % (string.replace(geo_Collection_ShortName,'-','_'))
+ geoXmlFile = path.join(ADL_HOME,'xml/VIIRS',geoXmlFile)
+ if path.exists(geoXmlFile):
+ LOG.debug("We are using for %s: %s,%s" %(geo_Collection_ShortName,geoXmlFile,geoFiles[0]))
+
+ # Open the geolocation blob and get the latitude and longitude
+
+ endian = self.sdrEndian
+
+ geoBlobObj = adl_blob.map(geoXmlFile,geoFiles[0], endian=endian)
+
+ # Get scan_mode to find any bad scans
+
+ scanMode = geoBlobObj.scan_mode[:]
+ badScanIdx = np.where(scanMode==254)[0]
+ LOG.debug("Bad Scans: %r" % (badScanIdx))
+
+ # Detemine the min, max and range of the latitude and longitude,
+ # taking care to exclude any fill values.
+
+ if longFormGeoNames :
+ if endian==adl_blob.BIG_ENDIAN:
+ latitude = getattr(geoBlobObj,'latitude').byteswap()
+ longitude = getattr(geoBlobObj,'longitude').byteswap()
+ latitude = latitude.astype('float')
+ longitude = longitude.astype('float')
+ else:
+ latitude = getattr(geoBlobObj,'latitude').astype('float')
+ longitude = getattr(geoBlobObj,'longitude').astype('float')
+ else :
+ latitude = getattr(geoBlobObj,'lat').astype('float')
+ longitude = getattr(geoBlobObj,'lon').astype('float')
+
+ latitude = ma.masked_less(latitude,-800.)
+ latMin,latMax = np.min(latitude),np.max(latitude)
+ latRange = latMax-latMin
+
+ longitude = ma.masked_less(longitude,-800.)
+ lonMin,lonMax = np.min(longitude),np.max(longitude)
+ lonRange = lonMax-lonMin
+
+ LOG.debug("min,max,range of latitide: %f %f %f" % (latMin,latMax,latRange))
+ LOG.debug("min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+
+ # Determine the latitude and longitude fill masks, so we can restore the
+ # fill values after we have scaled...
+
+ latMask = latitude.mask
+ lonMask = longitude.mask
+
+ # Check if the geolocation is in radians, convert to degrees
+ if 'RGEO' in geo_Collection_ShortName :
+ LOG.debug("Geolocation is in radians, convert to degrees...")
+ latitude = np.degrees(latitude)
+ longitude = np.degrees(longitude)
+
+ latMin,latMax = np.min(latitude),np.max(latitude)
+ latRange = latMax-latMin
+ lonMin,lonMax = np.min(longitude),np.max(longitude)
+ lonRange = lonMax-lonMin
+
+ LOG.debug("New min,max,range of latitude: %f %f %f" % (latMin,latMax,latRange))
+ LOG.debug("New min,max,range of longitude: %f %f %f" % (lonMin,lonMax,lonRange))
+
+ # Restore fill values to masked pixels in geolocation
+
+ geoFillValue = self.trimObj.sdrTypeFill['VDNE_FLOAT64_FILL'][latitude.dtype.name]
+ latitude = ma.array(latitude,mask=latMask,fill_value=geoFillValue)
+ self.latitude = latitude.filled()
+
+ geoFillValue = self.trimObj.sdrTypeFill['VDNE_FLOAT64_FILL'][longitude.dtype.name]
+ longitude = ma.array(longitude,mask=lonMask,fill_value=geoFillValue)
+ self.longitude = longitude.filled()
+
+ # Record the corners, taking care to exclude any bad scans...
+ nDetectors = 16
+ firstGoodScan = np.where(scanMode<=2)[0][0]
+ lastGoodScan = np.where(scanMode<=2)[0][-1]
+ firstGoodRow = firstGoodScan * nDetectors
+ lastGoodRow = lastGoodScan * nDetectors + nDetectors - 1
+
+ latCrnList = [latitude[firstGoodRow,0],latitude[firstGoodRow,-1],latitude[lastGoodRow,0],latitude[lastGoodRow,-1]]
+ lonCrnList = [longitude[firstGoodRow,0],longitude[firstGoodRow,-1],longitude[lastGoodRow,0],longitude[lastGoodRow,-1]]
+
+ # Check for dateline/pole crossings
+ num180Crossings = findDatelineCrossings(latCrnList,lonCrnList)
+ LOG.debug("We have %d dateline crossings."%(num180Crossings))
+
+ # Copy the geolocation information to the class object
+ self.latMin = latMin
+ self.latMax = latMax
+ self.latRange = latRange
+ self.lonMin = lonMin
+ self.lonMax = lonMax
+ self.lonRange = lonRange
+ self.scanMode = scanMode
+ self.latitude = latitude
+ self.longitude = longitude
+ self.latCrnList = latCrnList
+ self.lonCrnList = lonCrnList
+ self.num180Crossings = num180Crossings
+
+ # Parse the geolocation asc file to get struct information which will be
+ # written to the ancillary asc files
+
+ geoAscFileName = path.join(self.inDir,URID+".asc")
+ LOG.debug("\nOpening %s..." % (geoAscFileName))
+
+ geoAscFile = open(geoAscFileName,'rt')
+
+ self.RangeDateTimeStr = getAscLine(geoAscFile,"ObservedDateTime")
+ self.RangeDateTimeStr = string.replace(self.RangeDateTimeStr,"ObservedDateTime","RangeDateTime")
+ self.GRingLatitudeStr = getAscStructs(geoAscFile,"GRingLatitude",12)
+ self.GRingLongitudeStr = getAscStructs(geoAscFile,"GRingLongitude",12)
+
+ geoAscFile.close()
+
+
+ def _grid2Gran_bilinearInterp(self,dataLat, dataLon, gridData, gridLat, gridLon):
+ '''Granulates a gridded dataset using an input geolocation'''
+
+ nData = np.int64(dataLat.size)
+ gridRows = np.int32(gridLat.shape[0])
+ gridCols = np.int32(gridLat.shape[1])
+
+ data = np.ones(np.shape(dataLat),dtype=np.float64)* -999.9
+ dataIdx = np.ones(np.shape(dataLat),dtype=np.int64) * -99999
+
+ ANC_SCRIPTS_PATH = path.join(CSPP_RT_HOME,'viirs')
+
+ libFile = path.join(ANC_SCRIPTS_PATH,'libgriddingAndGranulation.so.1.0.1')
+ LOG.debug("Gridding and granulation library file: %s" % (libFile))
+ lib = ctypes.cdll.LoadLibrary(libFile)
+ grid2gran_bilinearInterp = lib.grid2gran_bilinearInterp
+ grid2gran_bilinearInterp.restype = None
+ grid2gran_bilinearInterp.argtypes = [
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ctypes.c_int64,
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_double,ndim=2,shape=(gridRows,gridCols),flags='C_CONTIGUOUS'),
+ ndpointer(ctypes.c_int64,ndim=1,shape=(nData),flags='C_CONTIGUOUS'),
+ ctypes.c_int32,
+ ctypes.c_int32
+ ]
+
+ '''
+ int snapGrid_ctypes(double *lat,
+ double *lon,
+ double *data,
+ long nData,
+ double *gridLat,
+ double *gridLon,
+ double *gridData,
+ long *gridDataIdx,
+ int nGridRows,
+ int nGridCols
+ )
+ '''
+
+ LOG.debug("Calling C routine grid2gran_bilinearInterp()...")
+
+ retVal = grid2gran_bilinearInterp(dataLat,
+ dataLon,
+ data,
+ nData,
+ gridLat,
+ gridLon,
+ gridData,
+ dataIdx,
+ gridRows,
+ gridCols)
+
+ LOG.debug("Returning from C routine grid2gran_bilinearInterp()")
+
+ return data,dataIdx
+
+
+ def granulate(self,ANC_objects):
+ '''
+ Granulate the ancillary dataset.
+ '''
+ LOG.info("Granulating %s ..." % (self.collectionShortName))
+
+ degInc = 0.5
+
+ lats = np.arange(361.)*degInc - 90.
+ lons = np.arange(720.)*degInc - 180.
+ latitude = self.latitude
+ longitude = self.longitude
+
+ gridData = self.gridData[:,:]
+
+ if self.num180Crossings != 2 :
+
+ #gridData = np.roll(gridData,360) # old
+ gridData = np.roll(gridData,360,axis=1) # new
+ gridLon,gridLat = np.meshgrid(lons,lats)
+
+ LOG.debug("start,end NCEP Grid Latitude values : %f,%f"%(gridLat[0,0],gridLat[-1,0]))
+ LOG.debug("start,end NCEP Grid Longitude values : %f,%f"%(gridLon[0,0],gridLon[0,-1]))
+
+ else :
+
+ negLonIdx = np.where(lons<0)
+ lons[negLonIdx] += 360.
+ lons = np.roll(lons,360)
+ gridLon,gridLat = np.meshgrid(lons,lats)
+
+ longitudeNegIdx = np.where(longitude < 0.)
+ longitude[longitudeNegIdx] += 360.
+
+ LOG.debug("start,end NCEP Grid Latitude values : %f,%f"%(gridLat[0,0],gridLat[-1,0]))
+ LOG.debug("start,end NCEP Grid Longitude values : %f,%f"%(gridLon[0,0],gridLon[0,-1]))
+
+ LOG.debug("min of gridData = %r"%(np.min(gridData)))
+ LOG.debug("max of gridData = %r"%(np.max(gridData)))
+
+ t1 = time()
+ data,dataIdx = self._grid2Gran_bilinearInterp(np.ravel(latitude),
+ np.ravel(longitude),
+ gridData.astype(np.float64),
+ gridLat,
+ gridLon)
+ t2 = time()
+ elapsedTime = t2-t1
+ LOG.info("Granulation took %f seconds for %d points" % (elapsedTime,latitude.size))
+
+ data = data.reshape(latitude.shape)
+ dataIdx = dataIdx.reshape(latitude.shape)
+
+ LOG.debug("Shape of granulated %s data is %s" % (self.collectionShortName,np.shape(data)))
+ LOG.debug("Shape of granulated %s dataIdx is %s" % (self.collectionShortName,np.shape(dataIdx)))
+
+ # Moderate resolution trim table arrays. These are
+ # bool arrays, and the trim pixels are set to True.
+ modTrimMask = self.trimObj.createModTrimArray(nscans=48,trimType=bool)
+
+ # Fill the required pixel trim rows in the granulated NCEP data with
+ # the ONBOARD_PT_FILL value for the correct data type
+
+ fillValue = self.trimObj.sdrTypeFill['ONBOARD_PT_FILL'][self.dataType]
+ data = ma.array(data,mask=modTrimMask,fill_value=fillValue)
+ self.data = data.filled()
+
+
+ def shipOutToFile(self):
+ ''' Pass the current class instance to this Utils method to generate
+ a blob/asc file pair from the input ancillary data object.'''
+
+ return shipOutToFile(self)
diff --git a/colourMapTools.py b/colourMapTools.py
deleted file mode 100644
index ec047fce8962a63a507e91de95107118ca1d7b24..0000000000000000000000000000000000000000
--- a/colourMapTools.py
+++ /dev/null
@@ -1,133 +0,0 @@
-#!/usr/bin/env python
-
-import numpy as np
-from matplotlib.colors import Colormap, normalize, LinearSegmentedColormap
-from types import IntType, FloatType, ListType
-from scipy import interpolate
-
-def cmap_discretize(cmap, N):
- """Return a discrete colormap from the continuous colormap cmap.
-
- cmap: colormap instance, eg. cm.jet.
- N: Number of colors.
-
- Example
- x = resize(arange(100), (5,100))
- djet = cmap_discretize(cm.jet, 5)
- imshow(x, cmap=djet)
- """
-
- cdict = cmap._segmentdata.copy()
- # N colors
- colors_i = np.linspace(0,1.,N)
- # N+1 indices
- indices = np.linspace(0,1.,N+1)
- for key in ('red','green','blue'):
- # Find the N colors
- D = np.array(cdict[key])
- I = interpolate.interp1d(D[:,0], D[:,1])
- #I = np.interp(indices,D[:,0], D[:,1])
- colors = I(colors_i)
- # Place these colors at the correct indices.
- A = np.zeros((N+1,3), float)
- A[:,0] = indices
- A[1:,1] = colors
- A[:-1,2] = colors
- # Create a tuple for the dictionary.
- L = []
- for l in A:
- L.append(tuple(l))
- cdict[key] = tuple(L)
-
- # Return colormap object.
- return LinearSegmentedColormap('colormap',cdict,1024)
-
-def cmap_reverse(my_cmap):
- my_newdict = {}
- for items in my_cmap._segmentdata :
- item_tuple = my_cmap._segmentdata[items]
- new_item_tuple = []
- item_tuple_revIter = reversed(item_tuple)
- for segments in item_tuple :
- dummyTuple = item_tuple_revIter.next()
- new_item_tuple.append(tuple((1.-dummyTuple[0],dummyTuple[2],dummyTuple[1])))
- my_newdict[items] = tuple(new_item_tuple)
-
- return LinearSegmentedColormap('my_newcolormap',my_newdict,256)
-
-####################################################################
-### Inserts black segment into colourmap up to splitVal, ###
-### compressing rest of colourmap towards high end. Colours ###
-### in colourmap are preserved. ###
-####################################################################
-
-def cmap_preTruncate(my_cmap,splitVal):
- my_newdict = {}
- for items in my_cmap._segmentdata :
- item_tuple = my_cmap._segmentdata[items]
- new_item_tuple = []
- new_item_tuple.append(tuple((0.0,0.0,0.0)))
- new_item_tuple.append(tuple((splitVal,0.0,item_tuple[0][2])))
-
- item_tuple_Iter = iter(item_tuple)
- item_tuple_Iter.next()
-
- while 1 :
- try :
- dummyTuple = item_tuple_Iter.next()
- newVal = splitVal + dummyTuple[0]*(1.0 - splitVal)
- new_item_tuple.append(tuple((newVal,dummyTuple[2],dummyTuple[1])))
- except StopIteration :
- break
-
- my_newdict[items] = tuple(new_item_tuple)
-
- return LinearSegmentedColormap('my_newcolormap',my_newdict,256)
-
-####################################################################
-### Replaces colourmap upto splitVal with black, and rescales ###
-### colourmap so that the transition occurs at startVal. ###
-### Colours in original colourmap up to splitval are lost. ###
-####################################################################
-
-def cmap_preTruncate2(my_cmap,startVal,splitVal):
- my_newdict = {}
- for items in my_cmap._segmentdata :
- item_tuple = my_cmap._segmentdata[items]
- colourTuple = []
- colourTupleNew = []
-
- # Create tuple containing tuples with colorbar coord > splitVal
- item_tuple_Iter = iter(item_tuple)
- while 1 :
- try :
- dummyTuple = item_tuple_Iter.next()
- if (dummyTuple[0] > splitVal) :
- colourTuple.append([dummyTuple[0],dummyTuple[1],dummyTuple[2]])
- except StopIteration :
- break
-
- # Rescale colorbar coord to range [ startVal-->1.0 ]
- oldRange = colourTuple[-1][0] - colourTuple[0][0]
- newRange = 1.0 - startVal
-
- for tup in np.arange(np.shape(colourTuple)[0]) :
- newCbarCoord = newRange*(colourTuple[tup][0]-colourTuple[0][0])/oldRange + startVal
- colourTupleNew.append([newCbarCoord,colourTuple[tup][1],colourTuple[tup][2]])
-
- colourTupleNew[0][1] = 0.0
-
- new_item_tuple = []
- new_item_tuple.append((0.0,0.0,0.0))
- colourTupleNew_iter = iter(colourTupleNew)
-
- while 1 :
- try :
- new_item_tuple.append(colourTupleNew_iter.next())
- except StopIteration :
- break
-
- my_newdict[items] = tuple(new_item_tuple)
-
- return LinearSegmentedColormap('my_newcolormap',my_newdict,256)
-
diff --git a/griddingAndGranulation.c b/griddingAndGranulation.c
index ee9cb43e1ccb79e5ba1374f3581e2ae202d8e340..ac1dd476919f18e8bd0a0b11d64765322d2d7bcf 100644
--- a/griddingAndGranulation.c
+++ b/griddingAndGranulation.c
@@ -1,8 +1,8 @@
-/* file_Date = $Date$
- * file_Revision = $Revision$
- * file_Author = $Author$
- * file_HeadURL = $HeadURL$
- * file_Id = $Id$
+/* file_Date = $Date: 2014-09-08 17:09:53 -0500 (Mon, 08 Sep 2014) $
+ * file_Revision = $Revision: 2226 $
+ * file_Author = $Author: geoffc $
+ * file_HeadURL = $HeadURL: https://svn.ssec.wisc.edu/repos/jpss_adl/trunk/scripts/edr/griddingAndGranulation.c $
+ * file_Id = $Id: griddingAndGranulation.c 2226 2014-09-08 22:09:53Z geoffc $
*
* ________________________
* griddingAndGranulation.c
@@ -385,7 +385,7 @@ int gran2grid_moments(double *lat,
/*
* _________
- * grid2gran
+ * grid2gran_nearest
* _________
*
* This routine takes the regular grid gridData[nGridRows*nGridCols],
@@ -395,8 +395,8 @@ int gran2grid_moments(double *lat,
*
* Inputs:
*
- * gridLat : Non-regular latitude grid (nGridRows*nGridCols)
- * gridLon : Non-regular longitude grid (nGridRows*nGridCols)
+ * gridLat : Regular latitude grid (nGridRows*nGridCols)
+ * gridLon : Regular longitude grid (nGridRows*nGridCols)
* gridData : Data array on grid defined by gridLat and gridLon
* grid (nGridRows*nGridCols)
* lat : Non-regular latitude grid (nData)
@@ -413,17 +413,17 @@ int gran2grid_moments(double *lat,
*
*
*/
-int grid2gran(double *lat,
- double *lon,
- double *data,
- long nData,
- double *gridLat,
- double *gridLon,
- double *gridData,
- long *gridDataIdx,
- int nGridRows,
- int nGridCols
- )
+int grid2gran_nearest(double *lat,
+ double *lon,
+ double *data,
+ long nData,
+ double *gridLat,
+ double *gridLon,
+ double *gridData,
+ long *gridDataIdx,
+ int nGridRows,
+ int nGridCols
+ )
{
long idx;
@@ -439,30 +439,30 @@ int grid2gran(double *lat,
int crnrPt,snapCrnrPt;
int rowInBounds,colInBounds;
- printf("Shape of data is (%ld)\n",nData);
- printf("Shape of gridData is (%d, %d)\n", nGridRows,nGridCols);
- printf("TRUE = %d\n", TRUE);
- printf("FALSE = %d\n", FALSE);
+ /*printf("Shape of data is (%ld)\n",nData);*/
+ /*printf("Shape of gridData is (%d, %d)\n", nGridRows,nGridCols);*/
+ /*printf("TRUE = %d\n", TRUE);*/
+ /*printf("FALSE = %d\n", FALSE);*/
int numShow = 10;
int i;
- printf("\nDisplaying lat,lon,data...\n");
- for (i=0;i<numShow;i++){
- printf("%6.1f %6.1f %6.1f\n",lat[i],lon[i],data[i]);
- }
+ /*printf("\nDisplaying lat,lon,data...\n");*/
+ /*for (i=0;i<numShow;i++){*/
+ /*printf("%6.1f %6.1f %6.1f\n",lat[i],lon[i],data[i]);*/
+ /*}*/
// Determine the lat and lon grid spacings
- printf("nGridRows = %d\n", nGridRows);
- printf("nGridCols = %d\n", nGridCols);
- printf("gridLat[nGridCols] = %f\n", gridLat[nGridCols]);
- printf("gridLat[0] = %f\n", gridLat[0]);
- printf("gridLon[1] = %f\n", gridLon[1]);
- printf("gridLon[0] = %f\n", gridLon[0]);
+ /*printf("nGridRows = %d\n", nGridRows);*/
+ /*printf("nGridCols = %d\n", nGridCols);*/
+ /*printf("gridLat[nGridCols] = %f\n", gridLat[nGridCols]);*/
+ /*printf("gridLat[0] = %f\n", gridLat[0]);*/
+ /*printf("gridLon[1] = %f\n", gridLon[1]);*/
+ /*printf("gridLon[0] = %f\n", gridLon[0]);*/
gridLatInc = fabs(gridLat[nGridCols]-gridLat[0]);
gridLonInc = fabs(gridLon[1]-gridLon[0]);
- printf("gridLatInc,gridLonInc = (%6.1f %6.1f)\n",gridLatInc,gridLonInc);
+ /*printf("gridLatInc,gridLonInc = (%8.3f %8.3f)\n",gridLatInc,gridLonInc);*/
// Loop through non-gridded data points, find matching gridpoint, and assign
// gridpoint data value to that non-gridded data point.
@@ -553,3 +553,341 @@ int grid2gran(double *lat,
return(0);
}
+
+
+/*
+ * _________
+ * grid2gran_weightedAvg
+ * _________
+ *
+ * This routine takes the regular grid gridData[nGridRows*nGridCols],
+ * defined by gridLat[nGridRows*nGridCols] and gridLon[nGridRows*nGridRows],
+ * and regrids to the data[nData] array, defined on the non-regular grid
+ * defined by lat[ndata] and lon[nData].
+ *
+ * Inputs:
+ *
+ * gridLat : Non-regular latitude grid (nGridRows*nGridCols)
+ * gridLon : Non-regular longitude grid (nGridRows*nGridCols)
+ * gridData : Data array on grid defined by gridLat and gridLon
+ * grid (nGridRows*nGridCols)
+ * lat : Non-regular latitude grid (nData)
+ * lon : Non-regular longitude grid (nData)
+ * nData : Length of non-regular grid arrays
+ * nGridRows : Number of grid rows of the regular grid
+ * nGridCols : Number of grid columns of the regular grid
+ *
+ * Outputs:
+ *
+ * data : Data array on non-regular grid defined by lat and lon (nData)
+ * gridDataIdx : Array of indices into gridData, defined by lat and lon
+ * non-regular grid (nData)
+ *
+ *
+ */
+int grid2gran_weightedAvg(double *lat,
+ double *lon,
+ double *data,
+ long nData,
+ double *gridLat,
+ double *gridLon,
+ double *gridData,
+ long *gridDataIdx,
+ int nGridRows,
+ int nGridCols
+ )
+{
+
+ long idx;
+ int latGridIdxLo, latGridIdxHi, lonGridIdxLo, lonGridIdxHi;
+
+ double latVal,lonVal,dataVal;
+ double gridLatInc,gridLonInc;
+
+ int gridLatPt,gridLonPt;
+ int gridLatPoints[4], gridLonPoints[4];
+ double weight, weights[4], weightSum, weightedSum, eps;
+ double wAvgData;
+ double minDist,dist,latDist,lonDist;
+ double gridLatVal,gridLonVal;
+ int crnrPt,snapCrnrPt;
+ int rowInBounds,colInBounds;
+
+ /*printf("Shape of data is (%ld)\n",nData);*/
+ /*printf("Shape of gridData is (%d, %d)\n", nGridRows,nGridCols);*/
+ /*printf("TRUE = %d\n", TRUE);*/
+ /*printf("FALSE = %d\n", FALSE);*/
+
+ int numShow = 10;
+ int i;
+ /*printf("\nDisplaying lat,lon,data...\n");*/
+ /*for (i=0;i<numShow;i++){*/
+ /*printf("%6.1f %6.1f %6.1f\n",lat[i],lon[i],data[i]);*/
+ /*}*/
+
+ // Determine the lat and lon grid spacings
+ /*printf("nGridRows = %d\n", nGridRows);*/
+ /*printf("nGridCols = %d\n", nGridCols);*/
+ /*printf("gridLat[nGridCols] = %f\n", gridLat[nGridCols]);*/
+ /*printf("gridLat[0] = %f\n", gridLat[0]);*/
+ /*printf("gridLon[1] = %f\n", gridLon[1]);*/
+ /*printf("gridLon[0] = %f\n", gridLon[0]);*/
+
+ gridLatInc = fabs(gridLat[nGridCols]-gridLat[0]);
+ gridLonInc = fabs(gridLon[1]-gridLon[0]);
+
+ /*printf("gridLatInc,gridLonInc = (%6.1f %6.1f)\n",gridLatInc,gridLonInc);*/
+
+ // Loop through non-gridded data points, find matching gridpoint, and assign
+ // gridpoint data value to that non-gridded data point.
+
+ eps = 1.0e-12;
+
+ for (idx=0;idx<nData;idx++){
+ /*for (idx=0;idx<numShow;idx++){*/
+ latVal = lat[idx];
+ lonVal = lon[idx];
+
+ // Determine lat/lon grid indices which bound the non-gridded point
+ latGridIdxLo = (int) floor((latVal-gridLat[0])/gridLatInc);
+ latGridIdxHi = latGridIdxLo + 1;
+ lonGridIdxLo = (int) floor((lonVal-gridLon[0])/gridLonInc);
+ lonGridIdxHi = lonGridIdxLo + 1;
+
+ rowInBounds = TRUE;
+ colInBounds = TRUE;
+
+ // If the grid indices bounding the non-gridded point are off the
+ // grid, mark this non-gridded point as out-of-bounds.
+ if ((latGridIdxLo<0) || (latGridIdxHi>=nGridRows)){
+ rowInBounds = FALSE;
+ /*printf("Row idx out of bounds...\n");*/
+ }
+ if ((lonGridIdxLo<0) || (lonGridIdxHi>=nGridCols)){
+ colInBounds = FALSE;
+ /*printf("Column idx out of bounds...\n");*/
+ }
+
+ if (rowInBounds==FALSE){
+ continue;
+ }else if (colInBounds==FALSE){
+ continue;
+ }else{
+ gridLatPoints[0] = latGridIdxLo;
+ gridLatPoints[1] = latGridIdxLo;
+ gridLatPoints[2] = latGridIdxHi;
+ gridLatPoints[3] = latGridIdxHi;
+
+ gridLonPoints[0] = lonGridIdxLo;
+ gridLonPoints[1] = lonGridIdxHi;
+ gridLonPoints[2] = lonGridIdxLo;
+ gridLonPoints[3] = lonGridIdxHi;
+
+ minDist = 1.e+30;
+ snapCrnrPt = 0;
+
+ weightSum = 0.;
+ weightedSum = 0.;
+ wAvgData = 0.;
+
+ // Loop through the corners bounding the non-gridded point,
+ // to calculate the inverse-distance weights.
+ for (crnrPt=0;crnrPt<4;crnrPt++){
+
+ gridLatPt = (int) gridLatPoints[crnrPt];
+ gridLonPt = (int) gridLonPoints[crnrPt];
+
+ // Get the lat and lon values at this corner
+ gridLatVal = gridLat[nGridCols * gridLatPt + gridLonPt];
+ gridLonVal = gridLon[nGridCols * gridLatPt + gridLonPt];
+
+ // The Pythagorean distance of this corner from the non-gridded
+ // data point...
+ latDist = latVal-gridLatVal;
+ lonDist = lonVal-gridLonVal;
+ dist = sqrt(latDist*latDist + lonDist*lonDist);
+
+ weight = 1./(dist + eps);
+
+ weightedSum += weight * gridData[nGridCols * gridLatPt + gridLonPt];
+
+ weightSum += weight;
+ }
+
+ // Assign the value of the grid point closest to the non-gridded data
+ // point to that non-gridded data point .
+ data[idx] = weightedSum / weightSum;
+
+ // Save the index of the closest grid data point to the non-gridded
+ // data point.
+ gridDataIdx[idx] = nGridCols * gridLatPt + gridLonPt;
+
+ }
+ }
+
+ return(0);
+}
+
+
+/*
+ * _________
+ * grid2gran_bilinearInterp
+ * _________
+ *
+ * This routine takes the regular grid gridData[nGridRows*nGridCols],
+ * defined by gridLat[nGridRows*nGridCols] and gridLon[nGridRows*nGridRows],
+ * and regrids to the data[nData] array, defined on the non-regular grid
+ * defined by lat[ndata] and lon[nData].
+ *
+ * Inputs:
+ *
+ * gridLat : Non-regular latitude grid (nGridRows*nGridCols)
+ * gridLon : Non-regular longitude grid (nGridRows*nGridCols)
+ * gridData : Data array on grid defined by gridLat and gridLon
+ * grid (nGridRows*nGridCols)
+ * lat : Non-regular latitude grid (nData)
+ * lon : Non-regular longitude grid (nData)
+ * nData : Length of non-regular grid arrays
+ * nGridRows : Number of grid rows of the regular grid
+ * nGridCols : Number of grid columns of the regular grid
+ *
+ * Outputs:
+ *
+ * data : Data array on non-regular grid defined by lat and lon (nData)
+ * gridDataIdx : Array of indices into gridData, defined by lat and lon
+ * non-regular grid (nData)
+ *
+ *
+ * The notation below is based on the description at
+ * http://en.wikipedia.org/wiki/Bilinear_interpolation
+ *
+ * | | |
+ * |Q_12 |R_2 |Q_22
+ * lat_hi ------o----------o----------o----- \ \
+ * | | | | |
+ * | | | |-- dlat_plus |
+ * | | | | |
+ * latVal --|--------- P ---------|----- / |-- dlat
+ * | | | \ |
+ * | | | | |
+ * | | | |-- dlat_plus |
+ * |Q_11 |R_1 |Q_21 | |
+ * lat_lo ------o----------o----------o----- / /
+ * | | |
+ * | | |
+ * | lonVal |
+ * lon_lo lon_hi
+ *
+ * \__________/\_________/
+ * | |
+ * dlon_minus dlon_plus
+ *
+ * \_____________________/
+ * |
+ * dlon
+ *
+ *
+ */
+int grid2gran_bilinearInterp(double *lat,
+ double *lon,
+ double *data,
+ long nData,
+ double *gridLat,
+ double *gridLon,
+ double *gridData,
+ long *gridDataIdx,
+ int nGridRows,
+ int nGridCols
+ )
+{
+
+ long idx;
+ int latGridIdxLo, latGridIdxHi, lonGridIdxLo, lonGridIdxHi;
+
+ double latVal,lonVal,dataVal;
+ double gridLatInc,gridLonInc;
+ double lat_lo, lat_hi, dlat, dlat_minus, dlat_plus;
+ double lon_lo, lon_hi, dlon, dlon_minus, dlon_plus;
+ double Q_11, Q_12, Q_21, Q_22;
+ double f_R1, f_R2, f_P;
+
+ int gridLatPt,gridLonPt;
+ double eps;
+ double gridLatVal,gridLonVal;
+ int rowInBounds,colInBounds;
+
+ int numShow = 10;
+ int i;
+
+ // Calculate the grid lat and lon increments
+ gridLatInc = fabs(gridLat[nGridCols]-gridLat[0]);
+ gridLonInc = fabs(gridLon[1]-gridLon[0]);
+
+ // Loop through non-gridded data points, find the bounding gridpoints,
+ // bilinearly interpolate.
+
+ eps = 1.0e-12;
+
+ for (idx=0;idx<nData;idx++){
+
+ latVal = lat[idx];
+ lonVal = lon[idx];
+
+ // Determine lat/lon grid indices which bound the non-gridded point
+ latGridIdxLo = (int) floor((latVal-gridLat[0])/gridLatInc);
+ latGridIdxHi = latGridIdxLo + 1;
+ lonGridIdxLo = (int) floor((lonVal-gridLon[0])/gridLonInc);
+ lonGridIdxHi = lonGridIdxLo + 1;
+
+ rowInBounds = TRUE;
+ colInBounds = TRUE;
+
+ // If the grid indices bounding the non-gridded point are off the
+ // grid, mark this non-gridded point as out-of-bounds.
+ if ((latGridIdxLo<0) || (latGridIdxHi>=nGridRows)){
+ rowInBounds = FALSE;
+ }
+ if ((lonGridIdxLo<0) || (lonGridIdxHi>=nGridCols)){
+ colInBounds = FALSE;
+ }
+
+ if (rowInBounds==FALSE){
+ continue;
+ }else if (colInBounds==FALSE){
+ continue;
+ }else{
+
+ lat_lo = gridLat[nGridCols * latGridIdxLo + lonGridIdxLo];
+ lat_hi = gridLat[nGridCols * latGridIdxHi + lonGridIdxLo];
+ dlat_minus = latVal - lat_lo;
+ dlat_plus = lat_hi - latVal;
+ dlat = lat_hi - lat_lo;
+
+ lon_lo = gridLon[nGridCols * latGridIdxLo + lonGridIdxLo];
+ lon_hi = gridLon[nGridCols * latGridIdxLo + lonGridIdxHi];
+ dlon_minus = lonVal - lon_lo;
+ dlon_plus = lon_hi - lonVal;
+ dlon = lon_hi - lon_lo;
+
+ Q_11 = gridData[nGridCols * latGridIdxLo + lonGridIdxLo];
+ Q_12 = gridData[nGridCols * latGridIdxHi + lonGridIdxLo];
+ Q_21 = gridData[nGridCols * latGridIdxLo + lonGridIdxHi];
+ Q_22 = gridData[nGridCols * latGridIdxHi + lonGridIdxHi];
+
+ f_R1 = (dlon_plus * Q_11 + dlon_minus * Q_21) / dlon;
+ f_R2 = (dlon_plus * Q_12 + dlon_minus * Q_22) / dlon;
+ f_P = (dlat_plus * f_R1 + dlat_minus * f_R2) / dlat;
+
+ // Assign the value of the grid point closest to the non-gridded data
+ // point to that non-gridded data point .
+ data[idx] = f_P;
+
+ // Save the index of the closest grid data point to the non-gridded
+ // data point.
+ gridDataIdx[idx] = nGridCols * gridLatPt + gridLonPt;
+
+ }
+ }
+
+ return(0);
+}
diff --git a/netCDF2bin_imager.py b/netCDF2bin_imager.py
deleted file mode 100644
index 2e0850563212421e3f16e81f03f4416973c42e9f..0000000000000000000000000000000000000000
--- a/netCDF2bin_imager.py
+++ /dev/null
@@ -1,568 +0,0 @@
-#!/usr/bin/env python
-# encoding: utf-8
-"""
-netCDF2bin_imager.py
-
-Purpose: Read band arrays from a netcdf4/cf created by pytroll/mpop, and write them to flat
-binary files.
-
-Input:
- * netCDF file containing band arrays
-
-Output:
- * flat binary files for each band
-
-Details:
- *
-
-Preconditions:
- * netCDF4 python module
-
-Optional:
- *
-
-Minimum commandline:
-
- python netCDF2bin_imager.py --input_files=INPUTFILES
-
-where...
-
- INPUTFILES: The fully qualified path to the input files. May be a directory or a file glob.
-
-
-Created by Geoff Cureton <geoff.cureton@ssec.wisc.edu> on 2014-03-25.
-Copyright (c) 2014 University of Wisconsin Regents. All rights reserved.
-
- This program is free software: you can redistribute it and/or modify
- it under the terms of the GNU General Public License as published by
- the Free Software Foundation, either version 3 of the License, or
- (at your option) any later version.
-
- This program is distributed in the hope that it will be useful,
- but WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- GNU General Public License for more details.
-
- You should have received a copy of the GNU General Public License
- along with this program. If not, see <http://www.gnu.org/licenses/>.
-"""
-
-file_Date = '$Date: 2014-01-23 09:02:22 -0600 (Thu, 23 Jan 2014) $'
-file_Revision = '$Revision: 1885 $'
-file_Author = '$Author: kathys $'
-file_HeadURL = '$HeadURL: https://svn.ssec.wisc.edu/repos/jpss_adl/trunk/scripts/edr/adl_viirs_edr.py $'
-file_Id = '$Id: adl_viirs_edr.py 1885 2014-01-23 15:02:22Z kathys $'
-
-__author__ = 'Geoff Cureton <geoff.cureton@ssec.wisc.edu>'
-__version__ = '$Id: adl_viirs_edr.py 1885 2014-01-23 15:02:22Z kathys $'
-__docformat__ = 'Epytext'
-
-import os, sys, logging, traceback
-from os import path,uname,environ
-import string
-import re
-import uuid
-from shutil import rmtree,copyfile
-from glob import glob
-from time import time
-from datetime import datetime,timedelta
-
-import numpy as np
-from numpy import ma
-import copy
-
-import matplotlib
-import matplotlib.cm as cm
-from matplotlib.colors import ListedColormap
-from matplotlib.figure import Figure
-
-matplotlib.use('Agg')
-from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
-
-# This must come *after* the backend is specified.
-import matplotlib.pyplot as ppl
-
-from mpl_toolkits.basemap import Basemap
-
-from netCDF4 import Dataset
-from netCDF4 import num2date
-
-
-def listData(file_attrs,file_dicts):
- '''
- List various file and data attributes
- '''
-
- print "\n### File Attributes ###\n"
- for keys in file_attrs.keys():
- print "{} : {}".format(keys,file_attrs[keys])
-
- print "\n### File Dimensions ###\n"
- print "Keys = {}".format(file_dicts['dimensions'].keys())
- for keys in file_dicts['dimensions'].keys():
- print "{} : {}\n".format(keys,file_dicts['dimensions'][keys])
-
- print "\n### File Variables ###\n"
- var_attrs = ['_FillValue','scale_factor','add_offset','coordinates','long_name','resolution','units','ndim','size','dtype','shape','dimensions']
- print "Keys = {}".format(file_dicts['variables'].keys())
- for var_key in file_dicts['variables'].keys():
- print "\nvar = {}".format(var_key)
- var = file_dicts['variables'][var_key]
- for attr in var_attrs:
- try:
- var_attr = getattr(var,attr)
- print "{}.{} = {}".format(var_key,attr,var_attr)
- except:
- pass
-
-
-
-def plotData(data,data_mask,units,pngName,stride=1,dpi=200):
- '''
- Plot the input dataset in swath projection
- '''
-
- # General Setup
- dataShape = data.shape
- aspectRatio = float(dataShape[0])/float(dataShape[1])
-
- figWidth,figHeight = 4.,aspectRatio*4.
- ax_rect = [0.05, 0.15, 0.90, 0.80 ] # [left,bottom,width,height]
-
- fig = Figure(figsize=(figWidth,figHeight))
- canvas = FigureCanvas(fig)
-
- ax = fig.add_axes(ax_rect)
-
- data = ma.array(data,mask=data_mask)
-
- im = ax.imshow(data,interpolation='nearest')
-
- if units == '%':
- txt = ax.set_title('Reflectance')
- elif units == 'K':
- txt = ax.set_title('Brightness Temperature')
- else:
- txt = ax.set_title('scaled data')
-
-
- ppl.setp(ax.get_xticklines(),visible=False)
- ppl.setp(ax.get_yticklines(),visible=False)
- ppl.setp(ax.get_xticklabels(), visible=False)
- ppl.setp(ax.get_yticklabels(), visible=False)
-
- ax.set_aspect('equal')
-
- cax_rect = [0.05 , 0.05, 0.90 , 0.05 ] # [left,bottom,width,height]
- cax = fig.add_axes(cax_rect,frameon=False) # setup colorbar axes
- cb = fig.colorbar(im, cax=cax, orientation='horizontal')
-
- if units == '%':
- txt = cax.set_title('Reflectance (%)')
- elif units == 'K':
- txt = cax.set_title('Brightness Temperature (K)')
- else:
- txt = cax.set_title('scaled units')
-
- # Redraw the figure
- canvas.draw()
- canvas.print_figure(pngName,dpi=dpi)
-
-
-def plotMapData(lat,lon,data,data_mask,units,pngName,stride=1,pointSize=1.,scatterPlot=False, \
- lat_0=None,lon_0=None,latMin=None,lonMin=None,latMax=None,lonMax=None,dpi=200):
- '''
- Plot the input dataset in mapped to particular projection
- '''
-
- # General Setup
- figWidth,figHeight = 8.,8.
- ax_rect = [0.05, 0.15, 0.90, 0.80 ] # [left,bottom,width,height]
-
- fig = Figure(figsize=(figWidth,figHeight))
- canvas = FigureCanvas(fig)
-
- ax = fig.add_axes(ax_rect)
-
- lat_0,lon_0 = 0.,0.
-
- #print "lonMin = {} degrees".format(lonMin)
- #print "lonMax = {} degrees".format(lonMax)
- #print "latMin = {} degrees".format(latMin)
- #print "latMax = {} degrees".format(latMax)
-
- m = Basemap(projection='cyl',lon_0=lon_0,lat_0=lat_0,lat_ts=lat_0,ax=ax,
- llcrnrlon=lonMin,
- llcrnrlat=latMin,
- urcrnrlon=lonMax,
- urcrnrlat=latMax
- )
-
- x,y=m(lon[::stride,::stride],lat[::stride,::stride])
-
- m.drawcoastlines()
- #m.drawstates()
- m.drawcountries()
- m.fillcontinents(color='0.85',zorder=0)
- m.drawparallels(np.arange( -90, 91,30), color = '0.25', linewidth = 0.5,labels=[1,0,1,0])
- m.drawmeridians(np.arange(-180,180,30), color = '0.25', linewidth = 0.5,labels=[1,0,1,0])
-
- data = ma.array(data[::stride,::stride],mask=data_mask[::stride,::stride])
-
- if scatterPlot:
- cs = m.scatter(x,y,s=pointSize,c=data,axes=ax,edgecolors='none')
- else:
- cs = m.pcolor(x,y,data,axes=ax,edgecolors='none',antialiased=False)
-
- if units == '%':
- txt = ax.set_title('Reflectance')
- elif units == 'K':
- txt = ax.set_title('Brightness Temperature')
- else:
- txt = ax.set_title('scaled data')
-
- ppl.setp(ax.get_xticklines(),visible=False)
- ppl.setp(ax.get_yticklines(),visible=False)
- ppl.setp(ax.get_xticklabels(), visible=False)
- ppl.setp(ax.get_yticklabels(), visible=False)
-
- #ax.set_aspect('equal')
-
- cax_rect = [0.05 , 0.05, 0.90 , 0.05 ] # [left,bottom,width,height]
- cax = fig.add_axes(cax_rect,frameon=False) # setup colorbar axes
- #cb = fig.colorbar(im, cax=cax, orientation='horizontal')
- cb = fig.colorbar(cs, cax=cax, orientation='horizontal')
-
- if units == '%':
- txt = cax.set_title('Reflectance (%)')
- elif units == 'K':
- txt = cax.set_title('Brightness Temperature (K)')
- else:
- txt = cax.set_title('scaled units')
-
- # Redraw the figure
- canvas.draw()
- canvas.print_figure(pngName,dpi=dpi)
-
-
-class dataset():
- '''
- Define a "dataset" datatpe, which contains useful attributes for a particular
- dataset.
- '''
- def __init__(self):
- self.bandNum
- self.units = None
- self.dpi = 200
- self.arrayName = None
- self.mapRes = None
-
-
-def _argparse():
- '''
- Method to encapsulate the option parsing and various setup tasks.
- '''
-
- import argparse
-
- endianChoices = ['little','big']
-
- dataChoices=['REF1','REF2','REF3','BT1','BT2','REF','BT','ALL']
-
- mapRes = ['c','l','i']
-
- defaults = {
- 'input_file':None,
- 'dataset':'ALL',
- 'stride':1,
- 'unscaled':False,
- 'scatter_plot':False,
- 'dpi':200,
- 'list_attrs':False,
- 'pointSize':1,
- 'make_plots':False,
- 'map_plot':False,
- 'no_binaries':False
- }
-
- description = '''Read band arrays from a netcdf4/cf created by pytroll/mpop, and write them to flat
-binary files.'''
- usage = "usage: %prog [mandatory args] [options]"
- version = __version__
-
- parser = argparse.ArgumentParser()
-
- # Mandatory arguments
-
- parser.add_argument('-i','--input_file',
- action='store',
- dest='input_file',
- type=str,
- required=True,
- help="The fully qualified path to the input file."
- )
-
- # Optional arguments
-
- parser.add_argument('--dset',
- action="store",
- dest="dataset",
- default=defaults["dataset"],
- type=str,
- choices=dataChoices,
- help='''The imager dataset(s) toVIIRS algorithm to run.\n\n
- Possible values are...
- {}.
- '''.format(dataChoices.__str__()[1:-1])
- )
-
- parser.add_argument('--unscaled',
- action="store_true",
- dest="doUnscale",
- default=defaults["unscaled"],
- help="Attempt to unscale the band arrays before output to flat binary files. [default: {}]".format(defaults["unscaled"])
- )
-
- parser.add_argument('--make_plots',
- action="store_true",
- dest="doPlots",
- default=defaults["make_plots"],
- help="Plot the band arrays to png files. [default: {}]".format(defaults["make_plots"])
- )
-
- parser.add_argument('--map_plot',
- action="store_true",
- dest="doMapPlots",
- default=defaults["map_plot"],
- help="Plot the band arrays to png files, using a map projection (has no effect if --make_plots is not set). [default: {}]".format(defaults["map_plot"])
- )
-
- parser.add_argument('--no_binaries',
- action="store_true",
- dest="noBinaries",
- default=defaults["no_binaries"],
- help="Plot the band arrays to png files. [default: {}]".format(defaults["make_plots"])
- )
-
- parser.add_argument('--plotMin',
- action="store",
- dest="plotMin",
- type=float,
- help="Minimum value to plot."
- )
-
- parser.add_argument('--plotMax',
- action="store",
- dest="plotMax",
- type=float,
- help="Maximum value to plot."
- )
-
- parser.add_argument('-d','--dpi',
- action="store",
- dest="dpi",
- default='200.',
- type=float,
- help="The resolution in dots per inch of the output png file. [default: {}]".format(defaults["dpi"])
- )
-
- parser.add_argument('-S','--stride',
- action="store",
- dest="stride",
- default=defaults["stride"],
- type=int,
- help="Sample every STRIDE pixels in the band data. [default: {}]".format(defaults["stride"])
- )
-
- parser.add_argument('--scatter_plot',
- action="store_true",
- dest="doScatterPlot",
- default=defaults["scatter_plot"],
- help="Generate the plot using a scatterplot approach."
- )
-
- parser.add_argument('--list_attrs',
- action="store_true",
- dest="doListAttrs",
- default=defaults["list_attrs"],
- help="List some file and dataset attributes."
- )
-
- parser.add_argument('-P','--pointSize',
- action="store",
- dest="pointSize",
- default=defaults["pointSize"],
- type=float,
- help="Size of the plot point used to represent each pixel. [default: {}]".format(defaults["pointSize"])
- )
-
- parser.add_argument('--lat_0',
- action="store",
- dest="lat_0",
- type=float,
- help="Center latitude of plot."
- )
-
- parser.add_argument('--lon_0',
- action="store",
- dest="lon_0",
- type=float,
- help="Center longitude of plot."
- )
-
- parser.add_argument('--latMin',
- action="store",
- dest="latMin",
- type=float,
- help="Minimum latitude to plot."
- )
-
- parser.add_argument('--latMax',
- action="store",
- dest="latMax",
- type=float,
- help="Maximum latitude to plot."
- )
-
- parser.add_argument('--lonMin',
- action="store",
- dest="lonMin",
- type=float,
- help="Minimum longitude to plot."
- )
-
- parser.add_argument('--lonMax',
- action="store",
- dest="lonMax",
- type=float,
- help="Maximum longitude to plot."
- )
-
- args = parser.parse_args()
-
- return args
-
-def main():
-
- options = _argparse()
-
- stride = options.stride
- lat_0 = options.lat_0
- lon_0 = options.lon_0
- latMin = options.latMin
- lonMin = options.lonMin
- latMax = options.latMax
- lonMax = options.lonMax
- doScatterPlot = options.doScatterPlot
- pointSize = options.pointSize
- dpi = options.dpi
-
- # Create the dataset dictionaries
- #dataDict = {}
- #if options.dataset == "REF":
-
- #if options.dataset == "BT":
- #if options.dataset == "ALL":
-
- dataDict = {
- 'REF1',
- 'REF2',
- 'REF3',
- 'BT1',
- 'BT2',
- 'REF',
- 'BT',
- 'ALL'
- }
-
- # Open the file object
- print "Opening {} ...".format(options.input_file)
- fileObj = Dataset(options.input_file)
-
- file_attrs = {}
- file_attrs['history'] = fileObj.history
- file_attrs['orbit'] = fileObj.orbit
- file_attrs['platform'] = fileObj.platform
-
- file_dicts = {}
- file_dicts['dimensions'] = fileObj.dimensions
- file_dicts['variables'] = fileObj.variables
-
- # List various file and data attributes
- if options.doListAttrs:
- listData(file_attrs,file_dicts)
-
- if not (options.noBinaries and not options.doPlots):
-
- print "We are reading the file Datasets"
-
- # Output the lats and lons
-
- lats = file_dicts['variables']['lat0'][:,:].astype(np.float32)
- lons = file_dicts['variables']['lon0'][:,:].astype(np.float32)
-
- if not options.noBinaries:
- geo_fileName = string.replace(options.input_file,'.nc','_latitude.bin')
- print "Writing latitude to {}...".format(geo_fileName)
- lats.tofile(geo_fileName)
-
- geo_fileName = string.replace(options.input_file,'.nc','_longitude.bin')
- print "Writing longitude to {}...".format(geo_fileName)
- lons.tofile(geo_fileName)
-
- # Output the bands
-
- for band_var in ['Image0','Image1']:
- #for band_var in ['Image0']:
- nBands = file_dicts['variables'][band_var].shape[2]
- dtype = file_dicts['variables'][band_var].dtype.name
- units = ''
-
- for band in range(nBands):
- #for band in [0]:
- var = file_dicts['variables'][band_var][:,:,band]
- fillValue = file_dicts['variables'][band_var]._FillValue
- print "fill value is {}".format(fillValue)
- if band_var == 'Image1':
- var_mask = ma.masked_inside(var,32720,32768).mask
- else:
- var_mask = ma.masked_equal(var,fillValue).mask
-
- scale_factor = file_dicts['variables'][band_var].scale_factor[band]
- offset = file_dicts['variables'][band_var].add_offset[band]
- print "{} band {} (scale,offset) = ({},{})".format(band_var,band,scale_factor,offset)
-
- if options.doUnscale :
- dtype = 'float32'
- units = file_dicts['variables'][band_var].units
- var = var * scale_factor + offset
- var = ma.array(var,mask=var_mask,fill_value=-999.9).filled()
-
- print "{} band {} fill value = {}".format(band_var,band,np.min(var))
- print "{} band {} min value = {}".format(band_var,band,np.min(ma.array(var,mask=var_mask)))
- print "{} band {} max value = {}".format(band_var,band,np.max(ma.array(var,mask=var_mask)))
-
- bin_fileName = string.replace(options.input_file,'.nc','')
- bin_fileName = "{}_{}_band{}_{}.bin".format(bin_fileName,band_var,band,dtype)
-
- if not options.noBinaries:
- print "Writing {} band {} to {}...".format(band_var,band,bin_fileName)
- var.astype(dtype).tofile(bin_fileName)
-
- if options.doPlots:
- png_fileName = string.replace(bin_fileName,'_{}.bin'.format(dtype),'.png')
- print "Plotting {} band {} to {}...".format(band_var,band,png_fileName)
- if options.doMapPlots:
- plotMapData(lats,lons,var,var_mask,units,png_fileName,stride=stride,lat_0=lat_0,lon_0=lon_0,latMin=latMin,lonMin=lonMin,latMax=latMax,lonMax=lonMax,scatterPlot=doScatterPlot,pointSize=pointSize,dpi=dpi)
- else:
- plotData(var,var_mask,units,png_fileName,stride=stride,dpi=dpi)
-
-
- print "Closing {} ...".format(options.input_file)
- fileObj.close()
-
- sys.exit(0)
-
-
-if __name__=='__main__':
- sys.exit(main())
diff --git a/resampleTools.py b/resampleTools.py
deleted file mode 100644
index cdeb81bc7654c43757bee9165a2fe3bc30f4de21..0000000000000000000000000000000000000000
--- a/resampleTools.py
+++ /dev/null
@@ -1,106 +0,0 @@
-import numpy as np
-import scipy.interpolate
-import scipy.ndimage
-
-def congrid(a, newdims, method='linear', centre=False, minusone=False):
- '''Arbitrary resampling of source array to new dimension sizes.
- Currently only supports maintaining the same number of dimensions.
- To use 1-D arrays, first promote them to shape (x,1).
-
- Uses the same parameters and creates the same co-ordinate lookup points
- as IDL''s congrid routine, which apparently originally came from a VAX/VMS
- routine of the same name.
-
- method:
- neighbour - closest value from original data
- nearest and linear - uses n x 1-D interpolations using
- scipy.interpolate.interp1d
- (see Numerical Recipes for validity of use of n 1-D interpolations)
- spline - uses ndimage.map_coordinates
-
- centre:
- True - interpolation points are at the centres of the bins
- False - points are at the front edge of the bin
-
- minusone:
- For example- inarray.shape = (i,j) & new dimensions = (x,y)
- False - inarray is resampled by factors of (i/x) * (j/y)
- True - inarray is resampled by(i-1)/(x-1) * (j-1)/(y-1)
- This prevents extrapolation one element beyond bounds of input array.
- '''
- if not a.dtype in [np.float64, np.float32]:
- a = np.cast[float](a)
-
- m1 = np.cast[int](minusone)
- ofs = np.cast[int](centre) * 0.5
- old = np.array( a.shape )
- ndims = len( a.shape )
- if len( newdims ) != ndims:
- print "[congrid] dimensions error. " \
- "This routine currently only support " \
- "rebinning to the same number of dimensions."
- return None
- newdims = np.asarray( newdims, dtype=float )
- dimlist = []
-
- if method == 'neighbour':
- for i in range( ndims ):
- base = np.indices(newdims)[i]
- dimlist.append( (old[i] - m1) / (newdims[i] - m1) \
- * (base + ofs) - ofs )
- cd = np.array( dimlist ).round().astype(int)
- newa = a[list( cd )]
- return newa
-
- elif method in ['nearest','linear']:
- # calculate new dims
- for i in range( ndims ):
- base = np.arange( newdims[i] )
- dimlist.append( (old[i] - m1) / (newdims[i] - m1) \
- * (base + ofs) - ofs )
- # specify old dims
- olddims = [np.arange(i, dtype = np.float) for i in list( a.shape )]
-
- # first interpolation - for ndims = any
- mint = scipy.interpolate.interp1d( olddims[-1], a, kind=method )
- newa = mint( dimlist[-1] )
-
- trorder = [ndims - 1] + range( ndims - 1 )
- for i in range( ndims - 2, -1, -1 ):
- newa = newa.transpose( trorder )
-
- mint = scipy.interpolate.interp1d( olddims[i], newa, kind=method )
- newa = mint( dimlist[i] )
-
- if ndims > 1 :
- # need one more transpose to return to original dimensions
- newa = newa.transpose( trorder )
-
- return newa
- elif method in ['spline']:
- oslices = [ slice(0,j) for j in old ]
- oldcoords = np.ogrid[oslices]
- nslices = [ slice(0,j) for j in list(newdims) ]
- newcoords = np.mgrid[nslices]
-
- newcoords_dims = range(np.rank(newcoords))
- #make first index last
- newcoords_dims.append(newcoords_dims.pop(0))
- newcoords_tr = newcoords.transpose(newcoords_dims)
- # makes a view that affects newcoords
-
- newcoords_tr += ofs
-
- deltas = (np.asarray(old) - m1) / (newdims - m1)
- newcoords_tr *= deltas
-
- newcoords_tr -= ofs
-
- newa = scipy.ndimage.map_coordinates(a, newcoords)
- return newa
- else :
- print "Congrid error: Unrecognized interpolation type.\n", \
- "Currently only \'neighbour\', \'nearest\',\'linear\',", \
- "and \'spline\' are supported."
- return None
-
diff --git a/snapgrid.m b/snapgrid.m
deleted file mode 100644
index b089158e3927b09853837d7f3b54193fd5187b6c..0000000000000000000000000000000000000000
--- a/snapgrid.m
+++ /dev/null
@@ -1,87 +0,0 @@
-function grid = snapgrid(spacing, geodata, param, first,grid)
-
-% input:
-% spacing - lat/lon grid spacing, e.g. 0.25, 0.5, 1.0, etc.
-% geodata.lat, geodata.lon, geodata.satzen, geodata.sunzen
-% param - the parameter you want to snap to grid, 'rad'for radiance or
-% 'cld' for cloud
-% first - 1=first; 0=not
-% grid - regrid out put from a previous run.
-% type - 1=day; 0=night
-% output:
-% regrid - data structure with regridded fields
-%
-% Nadia Smith, April 2010
-% Variation on Nick Bearson's code
-% =======================================================
-
-% if this is the first granule/file, then create new grids
-if first==1
- % day time
- grid.d_dist = ones(180/spacing, 360/spacing)*NaN;
- grid.d_satzen = ones(180/spacing, 360/spacing)*NaN;
- grid.d_solzen = ones(180/spacing, 360/spacing)*NaN;
- grid.d_param = ones(180/spacing, 360/spacing)*NaN;
- grid.d_type = ones(180/spacing, 360/spacing)*NaN;
- % night time
- grid.n_dist = ones(180/spacing, 360/spacing)*NaN;
- grid.n_satzen = ones(180/spacing, 360/spacing)*NaN;
- grid.n_solzen = ones(180/spacing, 360/spacing)*NaN;
- grid.n_param = ones(180/spacing, 360/spacing)*NaN;
- grid.n_type = ones(180/spacing, 360/spacing)*NaN;
-end
-
-for i=1:length(geodata.lat(:,1))
- for j=1:length(geodata.lat(1,:))
-
- % get the lat/lon of the data point
- lat = geodata.lat(i, j);
- lon = geodata.lon(i, j);
-
- if isnan(lat) == 0 && isnan(lon) == 0
-
- % adjust the lon/lat to be 0-360 and 0-180, and then
- % adjust them to our spacing grid for easier calculation
- adjlat = ((lat + 90) / spacing);
- adjlon = ((lon + 180) / spacing);
-
- % now find the closest grid crossing
- latbin = round(adjlat);
- lonbin = round(adjlon);
-
- % We don't want to reference the 0 spot, so we'll use the 180/360 instead (lat/lon)
- if latbin == 0
- latbin = 180 / spacing;
- end
- if lonbin == 0
- lonbin = 360 / spacing;
- end
-
- % determine our point's distance from the grid crossing
- % note that we're finding the hypotenuse distance, as that differs from just adding the two
- dist = sqrt(((latbin - adjlat)^2) + ((lonbin - adjlon)^2));
- if geodata.sunzen(i, j) < 84 % day time grid
- % now we're ready to get the current distance in that bin and see if ours is closer
- if (dist < (1/spacing)/4) && (isnan(grid.d_satzen(latbin, lonbin)) == 1 ||...
- abs(geodata.satzen(i, j)) < abs(grid.d_satzen(latbin, lonbin)))
- grid.d_dist(latbin, lonbin) = dist;
- grid.d_satzen(latbin, lonbin) = geodata.satzen(i, j);
- grid.d_sunzen(latbin, lonbin) = geodata.sunzen(i, j);
- grid.d_param(latbin, lonbin) = param(i, j);
- end
- else % night time grid
- % now we're ready to get the current distance in that bin and see if ours is closer
- if (dist < (1/spacing)/4) && (isnan(grid.n_satzen(latbin, lonbin)) == 1 ||...
- abs(geodata.satzen(i, j)) < abs(grid.n_satzen(latbin, lonbin)))
- grid.n_dist(latbin, lonbin) = dist;
- grid.n_satzen(latbin, lonbin) = geodata.satzen(i, j);
- grid.n_sunzen(latbin, lonbin) = geodata.sunzen(i, j);
- grid.n_param(latbin, lonbin) = param(i, j);
- end
- end
- end
- end
-end
-
-
-