[gdal-dev] Re: netcdf in gdal

[Norman Vine]

Frank Warmerdam writes:
> Peter Cromwell wrote:
> > Dear Frank,
> >
> > I am investigating the possibility of adding NetCDF to the
> > raster formats supported by GDAL.
> >
> > Are there any reasons why it would not be a good idea.
> >
> > To avoid duplicating effort, do you know whether someone has
> > done this already or is working on it?
> > >
> Peter,
>
> I think the folks who work on DODS have already written NetCDF support for
> GDAL. James, could you comment on that? What happened with regard to
> migrating your GDAL DODS driver back into the GDAL core?
>
> Also, others (ie. Norman Vine) have expressed an interest in having NetCDF
> support in GDAL.

Hi All

It would be nice to have native GDAL support for NetCDF
but it is relatively easy to use a Python bridge

Here is a Python snippet that will transform an arbritrary GDAL grid
into a NetCDF grid

Enjoy

Norman

==== cut ====
#!/usr/bin/env python
"""
GDAL_2_GMT -- Convert a GDAL file to a GMT grid file

built on top of Konrad Hinsen's NetCDF module
http://starship.python.net/crew/hinsen/scientific.html

Frank Warmerdam's GDAL extension module.
http://www.remotesensing.org/gdal/index.html

inspired the by wgrd function in Janko Hauser's NCfns module
http://starship.python.net/crew/jhauser/NCfns.py.html

Can optionally use NCAR's excellent regridding routines
http://esg.llnl.gov/cdat/
"""

__version__ = '0.2.1'
__email__ = 'nhv at c...'
__author__ = 'Norman Vine'

import os,sys,string
from Scientific.IO.NetCDF import *
import Numeric
import gdal
import math

# import adamsregrid

#---------------------
def wgrd(z,x,y,fname,title="no_title",name="dummy",
scale=1.,offset=0,
xunit="dummy",yunit="dummy",zunit="dummy"):
"""wgrd(z,x,y,fname,title=None,name=None,
scale=None,factor=None,
xunit=None,yunit=None,zunit=None)

Writes a grd-file suitable for GMT
"""

#print "NetCDFFile"
fn=NetCDFFile(fname,'w')

#print "Creating Fields"
# create all of our fields
fn.createDimension('side',2)
fn.createDimension('xysize',Numeric.product(z.shape))

fn.createVariable('x_range',Numeric.Float64,('side',))
fn.variables['x_range'].units=xunit

fn.createVariable('y_range',Numeric.Float64,('side',))
fn.variables['y_range'].units=yunit

fn.createVariable('z_range',Numeric.Float64,('side',))
fn.variables['z_range'].units=zunit

fn.createVariable('spacing',Numeric.Float64,('side',))
fn.createVariable('dimension',Numeric.Int32,('side',))

fn.createVariable('z',Numeric.Float32,('xysize',))
fn.variables['z'].long_name=name
fn.variables['z'].scale_factor=scale
fn.variables['z'].add_offset=offset
fn.variables['z'].node_offset=0

fn.title=title
fn.source="created by NumPy"

# Fill the fields
fn.variables['x_range'][0]=x[0]
fn.variables['x_range'][1]=x[1]
fn.variables['spacing'][0]=x[2]

fn.variables['y_range'][0]=y[0]
fn.variables['y_range'][1]=y[1]
fn.variables['spacing'][1]=y[2]

# determine zmin,zmax
print "checking min,max"
zmin = 100000.0
zmax = -100000.0
for zrow in z:
zmin = min(min(zrow),zmin)
zmax = max(max(zrow),zmax)
fn.variables['z_range'][0]=zmin
fn.variables['z_range'][1]=zmax

print "Writing array"
nsp=Numeric.asarray(z.shape)
fn.variables['dimension'][:]=nsp[::-1]
fn.variables['z'][:]=Numeric.ravel(z[:])
fn.close()
return


#---------------------
def gdal_2_gmt(filename):
"Convert GDAL grid to GMT grid."
z = []

# Open File with GDAL
ds = gdal.Open(filename)
if ds is None:
print 'Unable to open ', filename
sys.exit(1)

bands = ds.RasterCount
print "bands = %d"%(bands)
xsize = ds.RasterXSize
print "xsize = %d"%(xsize)
ysize = ds.RasterYSize
print "ysize = %d"%(ysize)
band_type = ds.GetRasterBand(1).DataType
print "band type ", band_type
# projection = ds.GetProjection()
# print "projection ",projection
geotransform = ds.GetGeoTransform()
print "geotransform ",geotransform
ulx = geotransform[0]
print "ulx %f"%(ulx)
uly = geotransform[3]
print "uly %f"%(uly)
lrx = ulx + geotransform[1] * xsize
print "lrx %f"%(lrx)
lry = uly + geotransform[5] * ysize
print "lry %f"%(lry)

xcellsize = geotransform[1]
print "xcellsize = %f"%(xcellsize)
ycellsize = -geotransform[5]
print "ycellsize = %f"%(ycellsize)

print "Reading Data"
z = ds.ReadAsArray()

# Test the NCAR's 'C' regrid extension :-))
# newxsize = 2049
# newysize = 2049
# print "newxsize = %d"%(newxsize)
# print "newysize = %d"%(newysize)
# r =
adamsregrid.Regrid(xsize,newxsize,'cubic',1,ysize,newysize,'cubic',0)
# z = r.rgrd(z)

# create Numeric array from Python array
# z = Numeric.array(z,Numeric.Float32)
# if z.shape[0] == 1 or z.shape[1] == 1:
# z = Numeric.ravel(z)

# xcellsize = math.fabs(lrx-ulx)/(newxsize-1)
# ycellsize = math.fabs(lry-uly)/(newysize-1)
# print "new xcellsize = %f"%(xcellsize)
# print "new ycellsize = %f"%(ycellsize)

x = [ ulx, lrx, xcellsize ]
y = [ lry, uly, ycellsize ]

# call our GMT grid file writer
basename,ext = os.path.splitext(filename)
filename = "%s.grd"%(basename)
print "Creating %s"%(filename)
wgrd(z,x,y,"%s"%(filename),title="no_title",name=basename,

scale=1.,offset=0,xunit="degrees",yunit="degrees",zunit="meters")

if __name__ == '__main__':
if len(sys.argv) == 2:
gdal_2_gmt(sys.argv[1])
else:
print "%s expecting 1 arg < filename to convert to GMT grid
>"%(sys.argv[0])