#! /usr/bin/env python # encoding: utf-8 ''' Created on August 14, 2011 Set up for _COLLECTION 6 FILES_ @author: nickb Copyright (c) 2011 University of Wisconsin SSEC. All rights reserved. For what the bits should be: http://modis-atmos.gsfc.nasa.gov/MOD35_L2/format.html ''' """ TEST: python ~/workspaces/python/pypre/src/convert_cloudout.py -i /Users/nickb/data/pypre/cloudmask_out_NEWTHRESH.hdf -o /Users/nickb/data/pypre/cloudmask_out_converted.hdf """ from files.HDF4File import HDF4File from files.HDF5File import HDF5File import os from optparse import OptionParser import numpy as np byte2bits = [''.join(['01'[i&(1<0] for b in xrange(7,-1,-1)]) for i in xrange(256)] def result_piecing(result, ouf): result0 = np.reshape(result[0], (2030, 1354)) result1 = np.reshape(result[1], (2030, 1354)) result2 = np.reshape(result[2], (2030, 1354)) result3 = np.reshape(result[3], (2030, 1354)) result4 = np.reshape(result[4], (2030, 1354)) result5 = np.reshape(result[5], (2030, 1354)) # BYTE 0 b0 = (result0 >> 0) & (2**0) b1 = (result0 >> 1) & (2**0 + 2**1) b3 = (result0 >> 3) & (2**0) b4 = (result0 >> 4) & (2**0) b5 = (result0 >> 5) & (2**0) b6 = (result0 >> 6) & (2**0 + 2**1) ouf.writeData("B00/Cloud Mask Flag", b0) ouf.writeData("B01/Unobstructed FOV Confidence Flag", b1) ouf.writeData("B03/DayNight Flag", b3) ouf.writeData("B04/Sunglint Flag", b4) ouf.writeData("B05/SnowIce Background Flag", b5) ouf.writeData("B06/LandWater Flag", b6) # BYTE 1 b8 = (result1 >> 0) & (2**0) b9 = (result1 >> 1) & (2**0) b10 = (result1 >> 2) & (2**0) b11 = (result1 >> 3) & (2**0) b12 = (result1 >> 4) & (2**0) b13 = (result1 >> 5) & (2**0) b14 = (result1 >> 6) & (2**0) b15 = (result1 >> 7) & (2**0) ouf.writeData("B08/Non-Cloud Obstruction Flag", b8 ) ouf.writeData("B09/Thin Cirrus Detected (Solar)", b9 ) ouf.writeData("B10/Snow Cover From Ancil Map", b10) ouf.writeData("B11/Thin Cirrus Detected (IR)", b11) ouf.writeData("B12/Cloud Adjacency", b12) ouf.writeData("B13/Cloud Flag (Ocean IR Threshold Test)", b13) ouf.writeData("B14/High Cloud Flag (CO2 Threshold est)", b14) ouf.writeData("B15/High Cloud Flag (6.7 um)", b15) # BYTE 2 b16 = (result2 >> 0) & (2**0) b17 = (result2 >> 1) & (2**0) b18 = (result2 >> 2) & (2**0) b19 = (result2 >> 3) & (2**0) b20 = (result2 >> 4) & (2**0) b21 = (result2 >> 5) & (2**0) b22 = (result2 >> 6) & (2**0) b23 = (result2 >> 7) & (2**0) ouf.writeData("B16/High Cloud Flag (1.38 um)", b16) ouf.writeData("B17/High Cloud Flag (3.7-12 um Test - night only)", b17) ouf.writeData("B18/Cloud Flag (IR Temperature Difference Tests)", b18) ouf.writeData("B19/Cloud Flag (3.9-11 um)", b19) ouf.writeData("B20/Cloud Flag (Visible Reflectance)", b20) ouf.writeData("B21/Cloud Flag (Visible Ratio)", b21) ouf.writeData("B22/Clearsky Restoral Test (NDVI in Coastal Areas)", b22) ouf.writeData("B23/Cloud Flag - Land and Polar Night (7.3-11 um)", b23) # BYTE 3 b24 = (result2 >> 0) & (2**0) b25 = (result2 >> 1) & (2**0) b26 = (result2 >> 2) & (2**0) b27 = (result2 >> 3) & (2**0) b28 = (result2 >> 4) & (2**0) b29 = (result2 >> 5) & (2**0) b30 = (result2 >> 6) & (2**0) b31 = (result2 >> 7) & (2**0) ouf.writeData("B24/Cloud Flag - Water 8.6-11 um Test", b24) ouf.writeData("B25/Clear-sky Restoral Test - Spatial Consistency (Ocean)", b25) ouf.writeData("B26/Clear-sky Restoral Test (polar night, land, sun-glint)", b26) ouf.writeData("B27/Cloud Flag - Surface Temperature Tests (water, night, land)", b27) ouf.writeData("B28/Suspended Dust Flag", b28) ouf.writeData("B29/Cloud Flag - Night Ocean 8.6 - 7.3 um Test", b29) ouf.writeData("B30/Cloud Flag - Night Ocean 11 um Variability Test", b30) ouf.writeData("B31/Cloud Flag - Night Ocean Low-Emissivity 3.9.11 um Test", b31) def cloudout_convert(infileName, navfileName, outfileName): print "Converting our cloud out file!" print "Output: ", outfileName inf = HDF5File(infileName) naf = HDF5File(navfileName) ouf = HDF5File(outfileName) latitude = naf.getDataset("latitude") longitude = naf.getDataset("longitude") latitude = np.reshape(latitude, (2030, 1354)) longitude = np.reshape(longitude, (2030, 1354)) ouf.writeData("Latitude", latitude) ouf.writeData("Longitude", longitude) result = inf.getDataset("result") result_piecing(result, ouf) def m35_convert(infileName, navfileName, outfileName): print "Converting M35 file!" print "Output: ", outfileName inf = HDF4File(infileName) # THE ONLY REAL DIFFERENCE HERE NOW naf = HDF5File(navfileName) ouf = HDF5File(outfileName) latitude = naf.getDataset("latitude") longitude = naf.getDataset("longitude") latitude = np.reshape(latitude, (2030, 1354)) longitude = np.reshape(longitude, (2030, 1354)) ouf.writeData("Latitude", latitude) ouf.writeData("Longitude", longitude) result = inf.getDataset("mod35/Data Fields/Swath Attributes/Cloud_Mask") result_piecing(result, ouf) def myinterp(data, ncols, nrows): from scipy import interpolate output = np.zeros((nrows, ncols)) rows, cols = data.shape xidxs = np.array(range(cols)) yidxs = np.array(range(rows)) z = data xx = np.linspace(xidxs.min(), xidxs.max(), ncols) yy = np.linspace(yidxs.min(), yidxs.max(), nrows) # Note that y and x are flipped? Weird, seems to work fine though. interpFunc = interpolate.RectBivariateSpline(yidxs, xidxs, z, kx=2, ky=2) output = interpFunc(yy, xx) return output 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.") parser.add_option("-n", "--navfile", action="store", type="string", dest="navfile", help="Navigation file.") parser.add_option("-m", action="store_true", dest="m35") parser.add_option("-c", action="store_false", dest="m35") 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 else: print "No input file provided." exit(1) if options.outfile is not None: outfile = options.outfile else: print "No output file provided." exit(1) if options.navfile is not None: navfile = options.navfile else: print "No nav file provided." exit(1) # Remove output file if it exists (so I can rapidly test!) if os.access(outfile, os.F_OK): os.remove(outfile) if options.m35 is None or options.m35 is False: cloudout_convert(infile, navfile, outfile) elif options.m35 is True: m35_convert(infile, navfile, outfile) else: print "Neither cloudmask output nor m35, WTF?" exit(1)