In [1]:

    

import os
from netCDF4 import Dataset
import numpy as np
from osgeo import gdal
from osgeo import osr
from pylab import rcParams
%matplotlib inline
import matplotlib.pyplot as plt


    

In [2]:

    

print "hello"


    

    




hello

In [7]:

    

grid = 'e2n_25'
#file = '/projects/PMESDR/vagrant/brodzik/' + grid + '/FD4a-E2T97-061-061.lis_dump.nc'
#file = '/projects/PMESDR/vagrant/brodzik/' + grid + '/FD4m-E2N97-061-061.lis_dump.nc'
file = '/projects/PMESDR/vagrant/brodzik/sirRSS/' + grid + '/FD4e-E2N97-061-061.lis_dump.nc'
f = Dataset( file, 'r', 'NETCDF4' )


    

In [9]:

    

# This should not be necessary, if the data are in the netCDF file properly
# Needs some investigation, could be FORTRAN-C-netCDF nonsense
tb = f.variables[ 'a_image' ][ : ]
cols, rows = shape( tb )
print cols, rows
tb = tb.reshape( rows, cols )


    

    




720 720

In [10]:

    

fig, ax = plt.subplots(1,1)
figsize(12,10)
ax.set_title( os.path.basename( file ) + ' (' + grid + ')' )
# Flip the image to display properly in python
plt.imshow(np.flipud(tb), cmap=plt.cm.gray, vmin=100, vmax=320 )
plt.axis('off')
#plt.colorbar( shrink=0.35, label='TB' )
plt.colorbar( shrink=0.95, label='TB' )


    

    
Out[10]:





<matplotlib.colorbar.Colorbar instance at 0x7f45d0e49a70>






    

In [11]:

    

outfile = file + '.new.nc'
driver = gdal.GetDriverByName('netCDF')


    

In [12]:

    

dst_ds = driver.Create(outfile, cols, rows, 1, gdal.GDT_UInt16)


    

In [13]:

    

proj = osr.SpatialReference()
# ease2_t:
# when we can connect to epsg v8.6, we should be able to do this: proj.SetWellKnownGeogCS("EPSG:6933")
#proj.SetFromUserInput("+proj=cea +lat_0=0 +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m")
# E2N
proj.SetFromUserInput("+proj=laea +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m")
# E2S
#proj.SetFromUserInput("+proj=laea +lat_0=-90 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m")
dst_ds.SetProjection(proj.ExportToWkt())


    

    
Out[13]:





0

In [14]:

    

print proj.ExportToPrettyWkt()


    

    




PROJCS["unnamed",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Lambert_Azimuthal_Equal_Area"],
    PARAMETER["latitude_of_center",90],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

In [15]:

    

# Thanks to web page at:
# http://geoexamples.blogspot.com/2012/01/creating-files-in-ogr-and-gdal-with.html
# The geotransform defines the relation between the raster coordinates x, y and the 
# geographic coordinates, using the following definition:
# Xgeo = geotransform[0] + Xpixel*geotransform[1] + Yline*geotransform[2]
# Ygeo = geotransform[3] + Xpixel*geotransform[4] + Yline*geotransform[5]
# The first and fourth parameters define the origin of the upper left pixel 
# The second and sixth parameters define the pixels size. 
# The third and fifth parameters define the rotation of the raster.
# These values are for EASE2_T25km only, values are meters
#E2T_25
#scale = 25025.26000
#E2T_3
#scale = 3128.15750
#E2N_25 or E2S_25
scale_x = 25000.00000
scale_y = -25000.00000
#E2N_3 or E2S_3
#scale = 3125.00000
# E2T
#map_UL_x = -17367530.44
#map_UL_y = -6756820.20000
# E2N or E2S
map_UL_x = -9000000.
map_UL_y = 9000000.
geotransform = (map_UL_x,scale_x,0.,map_UL_y,0.,scale_y)
dst_ds.SetGeoTransform(geotransform)


    

    
Out[15]:





0

In [16]:

    

# Convert tb array to UInt16
print tb.dtype
new = (tb + 0.5).astype(uint16)


    

    




float32

In [17]:

    

dst_ds.GetRasterBand(1).WriteArray(new)


    

    
Out[17]:





0

In [18]:

    

dst_ds = None


    

In [19]:

    

print outfile


    

    




/projects/PMESDR/vagrant/brodzik/sirRSS/e2n_25/FD4e-E2N97-061-061.lis_dump.nc.new.nc

In [20]:

    

print "Hello world"


    

    




Hello world

In [34]:

    




    

In [ ]: