Python + Fortran watermask resizing
'''
resize modis 44 w water mask using trs standard
scaling factors
usage:
resize_watermask.py <watermask_tile> <pixel resolution in metres>
example:
resize_watermask.py mod33h03v04.hdf 12
'''
import Nio
import numpy, pylab
import sys, os
import time
from arstats import ar_stats as mask
timestamp = time.strftime("%Y-%m-%d %H:%M:%S UTC", time.gmtime())
fpath = sys.argv[1]
dname = os.popen('dirname %s' % fpath).read().strip()
tile_no = os.popen('basename %s' % fpath).read().strip()
tile_no = tile_no.split(".")[2]
scale_f = sys.argv[2]
print ""
print "Output directory:\t%s" % dname
print "Tile:\t%s" % tile_no
odir = dname
oname = "mod44w.%s.%s.nc" % (tile_no, scale_f)
print "Output file path\t%s/%s" % (odir, oname)
print "Timestamp\t%s" % timestamp
scale = sys.argv[2] # pixel resolution
scf = float(scale)/250
iscf = int(scf)
# load water mask data
wm_file = Nio.open_file(fpath)
mask.wm = wm_file.variables["water_mask"][:]
ydim, xdim = mask.wm.shape
ndim = ydim / scf; mdim = xdim / scf
Y = int(ndim)
X = int(mdim)
mask.new_wm = numpy.ones((Y,X))*2
mask.scf = iscf
mask.mean()
# pylab.imshow(new_wm); pylab.show()
new_wm = mask.new_wm
new_wm = pylab.flipud(new_wm)
opath = "%s/%s" % (odir, oname)
if os.path.isfile(opath) == True: os.remove(opath)
ofile = Nio.open_file(opath, "c")
ofile.Conventions = "CF-1.4"
ofile.institute = "Max-Planck-Institute for Meteorology"
ofile.author = "Mikhail Itkin, mailto: mikhail.itkin@zmaw.de"
ofile.create_dimension("Y", Y) # latitude
ofile.create_dimension("X", X) # longitude
ofile.create_variable("water_mask", 'h', ("Y", "X"))
var = ofile.variables["water_mask"]
var.long_name = "Fraction of water area"
var.units = "percent"
var.source = "SRTM 0.9 km + MODIS 250m imagery"
var.references = "https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/land_water_mask_derived/land_water_mask_derived/mod44w"
ofile.comments = "%s $PROJECTS/watermask/resize_mask.py %s" % (timestamp, scf)
var._FillValue = numpy.array(-999,dtype=numpy.int16)
new_wm = numpy.where(new_wm > 1, -9.99, new_wm)
new_wm = numpy.array(new_wm*100, dtype = numpy.int16)
# new_wm = numpy.array(new_wm, dtype=numpy.uint16)
var[:] = new_wm
ofile.close()