#! /usr/bin/env python # encoding: utf-8 ''' Created on Jun 2, 2011 @author: nickb Copyright (c) 2011 University of Wisconsin SSEC. All rights reserved. ''' # FIXME: Messy PYTHONPATH so we can import from the directory above - is there a better way? import sys sys.path.append("..") from optparse import OptionParser #import glob import os #from os import path import numpy as np import ctypes from numpy.ctypeslib import ndpointer from variables.Variable import Variable from files.HDF5File import HDF5File #import _C_cloud_wrapper def run_cloudmask(lat, lon, senz, solz, ra, se, lsm, igbp, sst, nise, ndvi, eco, sfctemp, sfct_flag, tpw, modis, rowlen): # load the C extension library # libDir = path.dirname(__file__) # libDir = os.getcwd() libDir = "/data3/nickb/pypre/src/cloudmask" libFile = 'libcloudmask.so' libFile = "%s/%s" % (libDir, libFile) # print libFile # print "*****************" lib = ctypes.cdll.LoadLibrary(libFile) # print lib.__dict__ # print "*****************" # specify the ctypes argument types cloudwrapper_ctypes = lib.cm_wrapper cloudwrapper_ctypes.restype = ctypes.c_int32 # print "CREATING RESULT:" # print lat.shape # result = np.zeros((len(lat), 6), dtype=np.dtype("a1")); result = np.ones((6, (len(lat))), dtype=np.int8); qa = np.ones((6, (len(lat))), dtype=np.int8); cloudwrapper_ctypes.argtypes = [ ndpointer(ctypes.c_float,ndim=1,shape=lat.shape, flags='C_CONTIGUOUS'), #dlat ndpointer(ctypes.c_float,ndim=1,shape=lon.shape, flags='C_CONTIGUOUS'), #dlon ndpointer(ctypes.c_float,ndim=1,shape=senz.shape, flags='C_CONTIGUOUS'), #dsenZ ndpointer(ctypes.c_float,ndim=1,shape=solz.shape, flags='C_CONTIGUOUS'), #dsolZ ndpointer(ctypes.c_float,ndim=1,shape=ra.shape, flags='C_CONTIGUOUS'), #dra ndpointer(ctypes.c_float,ndim=1,shape=se.shape, flags='C_CONTIGUOUS'), #dse ndpointer(ctypes.c_float,ndim=1,shape=igbp.shape, flags='C_CONTIGUOUS'), #digbp ndpointer(ctypes.c_float,ndim=1,shape=sst.shape, flags='C_CONTIGUOUS'), #dsst ndpointer(ctypes.c_float,ndim=1,shape=ndvi.shape, flags='C_CONTIGUOUS'), #dndvi ndpointer(ctypes.c_float,ndim=1,shape=sfctemp.shape, flags='C_CONTIGUOUS'), #sfctemp ndpointer(ctypes.c_int32, ndim=1,shape=sfct_flag.shape, flags='C_CONTIGUOUS'), #sfct_flag ndpointer(ctypes.c_float,ndim=1,shape=tpw.shape, flags='C_CONTIGUOUS'), #tpw ndpointer(ctypes.c_float,ndim=2,shape=modis.shape, flags='C_CONTIGUOUS'), #dmodis ndpointer(ctypes.c_ubyte,ndim=1,shape=lsm.shape, flags='C_CONTIGUOUS'), #dlsm ndpointer(ctypes.c_ubyte,ndim=1,shape=nise.shape, flags='C_CONTIGUOUS'), #dnise ndpointer(ctypes.c_ubyte,ndim=1,shape=eco.shape, flags='C_CONTIGUOUS'), #deco ctypes.c_int32, ndpointer(ctypes.c_int8,ndim=2,shape=result.shape, flags='C_CONTIGUOUS'), ndpointer(ctypes.c_int8,ndim=2,shape=qa.shape, flags='C_CONTIGUOUS'), ] ''' The C function we're calling: unsigned char **cm_wrapper( float *dlat, float *dlon, float *dsenZ, float *dsolZ, float *dra, float *dse, float *digbp, float *dsst, float *dndvi, float *sfctemp, float *sfct_flag, float *tpw, float *dmodis, unsigned char *dlsm, unsigned char *dnise, unsigned char *deco, int lcol ) ''' retVal = cloudwrapper_ctypes(lat, lon, senz, solz, ra, se, igbp, sst, ndvi, sfctemp, sfct_flag, tpw, modis, lsm, nise, eco, rowlen, result, qa) return result, qa def main(ancilFile, outFile): h5ancil = HDF5File(ancilFile) latitude = Variable(sourceFile=h5ancil, varName='Latitude') latitude.getVariable() print "SHAPE:" print latitude.data.shape latitude.data = latitude.data.flatten().astype(np.float32) longitude = Variable(sourceFile=h5ancil, varName='Longitude') longitude.getVariable() longitude.data = longitude.data.flatten().astype(np.float32) SensorZenith = Variable(sourceFile=h5ancil, varName="SensorZenith") SensorZenith.getVariable() SensorZenith.data = SensorZenith.data.flatten().astype(np.float32) SolarZenith = Variable(sourceFile=h5ancil, varName='SolarZenith') SolarZenith.getVariable() SolarZenith.data = SolarZenith.data.flatten().astype(np.float32) RelativeAzimuth = Variable(sourceFile=h5ancil, varName='RelativeAzimuth') RelativeAzimuth.getVariable() RelativeAzimuth.data = RelativeAzimuth.data.flatten().astype(np.float32) SurfaceElevation = Variable(sourceFile=h5ancil, varName='SurfaceElevation') SurfaceElevation.getVariable() SurfaceElevation.data = SurfaceElevation.data.flatten().astype(np.float32) IGBPEco = Variable(sourceFile=h5ancil, varName='IGBP_LandCoverType') IGBPEco.getVariable() IGBPEco.data = IGBPEco.data.flatten().astype(np.float32) SST = Variable(sourceFile=h5ancil, varName='SeaSurfaceTemp') SST.getVariable() SST.data = SST.data.flatten().astype(np.float32) NDVI = Variable(sourceFile=h5ancil, varName='NDVI') NDVI.getVariable() NDVI.data = NDVI.data.flatten().astype(np.float32) sfctemp = Variable(sourceFile=h5ancil, varName='Surface Temperature sfct_flag') sfctemp.getVariable() sfctemp.data = sfctemp.data.flatten().astype(np.float32) tpw = Variable(sourceFile=h5ancil, varName='PrecipitableWater') tpw.getVariable() tpw.data = tpw.data / 10.0 # Divide by 10 as in alg29.cpp tpw.data = tpw.data.flatten().astype(np.float32) sfct_flag = Variable(sourceFile=h5ancil, varName='sfct_flag') sfct_flag.getVariable() sfct_flag.data = sfct_flag.data.flatten().astype(np.int32) # Char arrays get a different astype LandSeaMask = Variable(sourceFile=h5ancil, varName='MYD03_LandSeaMask') LandSeaMask.getVariable() LandSeaMask.data = LandSeaMask.data.flatten().astype(np.uint8) NISE = Variable(sourceFile=h5ancil, varName='NISE') NISE.getVariable() NISE.data = NISE.data.flatten().astype(np.uint8) ECO = Variable(sourceFile=h5ancil, varName="OlsonEcosystem") ECO.getVariable() ECO.data = ECO.data.flatten().astype(np.uint8) bandsMODISref = [1, 2, 3, 4, 5, 6, 7, 8, 9, 17, 18, 19, 26] modisData = np.zeros((36, len(latitude.data))) # modisData = [None] * 35 for band in bandsMODISref: modisVar = Variable(sourceFile=h5ancil, varName=('MODIS_reflectance_' + str(band))) modisVar.getVariable() # Adding a zeros array to make indexing simpler for us in C. # We are wasting memory by doing this, though - must find out # how much space "None" takes when it replaces a whole row or column. if modisVar.data is None: modisVar.data = np.zeros(latitude.data.size) modisData[band-1] = modisVar.data.flatten().astype(np.float32) bandsMODISbt = [20, 21, 22, 27, 28, 29, 31, 32, 33, 35] for band in bandsMODISbt: modisVar = Variable(sourceFile=h5ancil, varName=('MODIS_bt_' + str(band))) modisVar.getVariable() if modisVar.data is None: modisVar.data = np.zeros(latitude.data.size) modisData[band-1] = modisVar.data.flatten().astype(np.float32) # print modisData modisData = np.float32(modisData) # print modisData rowlen = len(modisData[0]) print "rowlen: ", rowlen print "********************************" print "REPRESENTATIVE MODISDATA VALUES!" for band in modisData: print band[0] result, qa = run_cloudmask(latitude.data, longitude.data, SensorZenith.data, SolarZenith.data, RelativeAzimuth.data, SurfaceElevation.data, LandSeaMask.data, IGBPEco.data, SST.data, NISE.data, NDVI.data, ECO.data, sfctemp.data, sfct_flag.data, tpw.data, modisData, rowlen) h5out = HDF5File(outFile) latitude = Variable(data=latitude.data, varName='latitude') h5out.writeVariable(latitude) longitude = Variable(data=longitude.data, varName='longitude') h5out.writeVariable(longitude) resultVar = Variable(data=result, varName='result') h5out.writeVariable(resultVar) qaVar = Variable(data=qa, varName='qa') h5out.writeVariable(qa) def handle_args(): parser = OptionParser() parser.add_option("-i", "--infile", action="store", type="string", dest="infile", help="Input file.") parser.add_option("-o", "--outFile", action="store", type="string", dest="outfile", help="Output file.") return parser.parse_args() if __name__ == '__main__': print "\n\n\n" # Handle arguments (options, args) = handle_args() if options.infile is not None: infile = options.infile print infile else: print "No input file provided." exit(1) if options.outfile is not None: outfile = options.outfile else: outfile = "cloudmask_out.hdf" if os.path.exists(outfile): os.remove(outfile) main(infile, outfile)