# Generate a GIS raster from the 8 NMQ data tiles! :)

import numpy
import mx.DateTime
import netCDF3
from PIL import Image
import os
import sys

# NW Corner of tiles
tiles = {
         1: [-130., 55.],
         2: [-110., 55.],
         3: [-90., 55.],
         4: [-80., 55.],
         5: [-130., 40.],
         6: [-110., 40.],
         7: [-90., 40.],
         8: [-80., 40.]
         }

def make_colorramp():
    """
    Make me a crude color ramp
    """
    c = numpy.zeros((256,3), numpy.int)
    
    # Gray for missing
    c[0,:] = [144,144,144]
    # Black to remove, eventually
    c[1,:] = [0,0,0]
    i = 2 
    for line in open('iem_lut256.tbl'):
      c[i,:] = map(int, line.split())
      i+=1
    return tuple( c.ravel() )

def make_fp(tile, gts):
    """
    Return a string for the filename expected for this timestamp
    """
    return "/mnt/a4/data/%s/nmq/tile%s/data/QPESUMS/grid/mosaic2d_nc/%s00.nc" % (
        gts.strftime("%Y/%m/%d"), tile, 
        gts.strftime("%Y%m%d-%H%M") )
    
def doit(gts):
    """
    Actually generate a PNG file from the 8 NMQ tiles
    """
    szx = 7000
    szy = 3500
    west = -130.
    north = 55.
    # Create the image data
    imgdata = numpy.zeros( (szy, szx), 'u1')
    # Loop over tiles
    for tile in range(1,9):
        fp = make_fp(tile, gts)
        if not os.path.isfile(fp):
            print "Missing", fp
            continue
        nc = netCDF3.Dataset( fp )
        val = nc.variables["hsr"][:,:] / 10.0 # convert to dBZ
        # -990 is no return
        # -9990 is missing
        # Bump up by one, so that we can set missing to color index 0
        #val += 1.0
        # Base is -30dBZ, so we add 60 to get us above zero, I hope
        #val = (val + 30.0) * 2.0
        val = numpy.where( val > 0., (val +30) * 2.0, val)
        val = numpy.where( val < -9900., 0., val)
        val = numpy.where( val < 0., 1., val)

        ysz, xsz = numpy.shape(val)
        x0 = (tiles[tile][0] - west) * 100.0
        y0 = (north - tiles[tile][1]) * 100.0
        #print tile, x0, xsz, y0, ysz, val[0,0]
        imgdata[y0:(y0+ysz-1),x0:(x0+xsz-1)] = val.astype('int')[:-1,:-1]
        nc.close()
    # Stress our color ramp
    #for i in range(256):
    #    imgdata[i*10:i*10+10,0:100] = i
    # Create Image
    png = Image.fromarray( imgdata )
    png.putpalette( make_colorramp() )
    png.save('q2.png')
    # Now we need to generate the world file
    o = open('/tmp/q2.wld', 'w')
    o.write("""   0.0100000000000%s
   0.00000
   0.00000
  -0.010000000000000000
%s
  %s""" % (gts.strftime("%Y%m%d%H%M"), west, north))
    o.close()
    # Inject WLD file
    pqstr = "/home/ldm/bin/pqinsert -p 'plot a %s bogus GIS/q2/hsr_%s.wld wld' /tmp/q2.wld" % (
                    gts.strftime("%Y%m%d%H%M"),gts.strftime("%Y%m%d%H%M") )
    os.system(pqstr)
    # Now we inject into LDM
    pqstr = "/home/ldm/bin/pqinsert -p 'plot ac %s gis/images/4326/q2/hsr.png GIS/q2/hsr_%s.png png' q2.png" % (
                    gts.strftime("%Y%m%d%H%M"),gts.strftime("%Y%m%d%H%M") )
    os.system(pqstr)
    # Create 900913 image
    cmd = "/mesonet/local/bin/gdalwarp -s_srs EPSG:4326 -t_srs EPSG:900913 -q -of GTiff -tr 1000.0 1000.0 q2.png q2.tif"
    os.system( cmd )
    # Insert into LDM
    pqstr = "/home/ldm/bin/pqinsert -p 'plot c %s gis/images/900913/q2/hsr.tif GIS/q2/hsr_%s.tif tif' q2.tif" % (
                    gts.strftime("%Y%m%d%H%M"),gts.strftime("%Y%m%d%H%M") )
    os.system(pqstr)
    
    os.unlink('q2.tif')

if __name__ == "__main__":
    if len(sys.argv) == 6:
        doit(mx.DateTime.DateTime(int(sys.argv[1]), int(sys.argv[2]), int(sys.argv[3]),
                                   int(sys.argv[4]), int(sys.argv[5])))
    else:
        gts = mx.DateTime.gmtime() - mx.DateTime.RelativeDateTime(minutes=10)
        offset = gts.minute % 5
        gts -= mx.DateTime.RelativeDateTime(minutes=offset)
        doit( gts )

@@ -1787,16 +1787,17 @@
 (val +30
+.
 ) * 2.0,
@@ -1839,17 +1839,16 @@
 l %3C -990
-0
 ., 0., v

divide by 10