#!/usr/bin/env python

"""Converts a land/water-mask in NetCDF (as output by grdlandmask) to image.

The image output is a 8-bit greyscale image with 0=land and 255=water.
Subsequent processing steps (e.g. with blur_lsm) can be used to create
"coastal" values.

Usage: lwm_netcdf2png levels.grd land_water_mask.png land_water_mask_latlon.nc

Example:

Depends on:
- ScientificPython to read/write NetCDF
- PIL to write imagefile
"""

# $Id$

import sys
import datetime
now = datetime.datetime.now
import Scientific.IO.NetCDF
import contextlib
import PIL.Image
##import matplotlib.pyplot
import numpy

class LandWaterMask(object):
    """Represents a land-water mask

    Attributes:

        lat  1-D ndarray with latitude values
        lon  1-D ndarray with longitude values
        lwm  2-D ndarray with 0=water, 1=land

    Methods:

        .__init__(f)
        .write(lat, lon, lwm)

    """

    def __init__(self, f):
        """Get land-water-mask from NetCDF (as output by grdlandmask).

        IN 

            f    path to NetCDF-file as output by grdlandmask

        OUT

            LandWaterMask object with attributes:
        """

        print now(), "Reading NetCDF"
        with contextlib.closing(Scientific.IO.NetCDF.NetCDFFile(f)) as nc:
            self.lon = nc.variables['x'][:]
            self.lat = nc.variables['y'][:]
            self.lwm = numpy.uint8(nc.variables['z'][:])

    def write(self, imfile, ncfile):
        """Write land-water-mask as image, coords in NetCDF

        IN

            imfile  path to image-file for land-water-mask
            ncfile  path to netcdf-file for coordinates

        """
        # write lsm
        print now(), "Writing image"
        imdata = self.lwm
        # change 1 to 255 for visibility
        imdata[imdata==1] = 255
        im = PIL.Image.fromarray(imdata)
        im.save(imfile)
        # write lat/lon
        print now(), "Writing lat/lon"
        with contextlib.closing(Scientific.IO.NetCDF.NetCDFFile(ncfile, 'w')) as nc:
            nc.createDimension("lat", self.lat.shape[0])
            nc.createDimension("lon", self.lon.shape[0])
            latv = nc.createVariable("lat", 'd', ("lat",))
            lonv = nc.createVariable("lon", 'd', ("lon",))
            latv[:] = self.lat
            lonv[:] = self.lon
            nc.flush()

def main():
    try:
        ncfile_in, imfile, ncfile_out = sys.argv[1:4]
    except ValueError:
        print >>sys.stderr, __doc__
        sys.exit(1)
    lwm = LandWaterMask(ncfile_in)
    lwm.write(imfile, ncfile_out)

if __name__ == "__main__":
    main()