Python + Fortran watermask resizing

{{{


#!/usr/bin/env python
'''
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()

}}}