#
# netCDF4 test script
#
# Rich Signell, original coding
# Curtis Price, small tweaks and doc
#
# Christoph Gohlke netCDF4 python library specially compiled
# for ArcGIS 10 at:
#   http://www.lfd.uci.edu/~gohlke/pythonlibs/
#     netCDF4-0.9.7-ArcGIS-10.0.win32-py2.6.exe
#     netCDF4-0.9.7-ArcGIS-10.1.win32-py2.7.exe

import os
import netCDF4
import numpy as np
import arcpy

# input netCDF url
url='http://geoport.whoi.edu/thredds/dodsC/bathy/crm_vol1.nc'
box = [-71.2,41,-70.2,42]
# output
arcpy.workspace  = "C:\\workspace"
rasterName = "traster"

outRaster = os.path.normpath(os.path.join(arcpy.workspace,rasterName))
nc=netCDF4.Dataset(url)
print "Source name: %s" % nc.title
lon=nc.variables['lon'][:]
lat=nc.variables['lat'][:]
bi=(lon>=box[0])&(lon<=box[2])
bj=(lat>=box[1])&(lat<=box[3])
zz=nc.variables['topo'][bj,bi]
xyCell = 3 / 3600. # 3-arc-second data
xyOrig = arcpy.Point(-71.2,41)
grid1=arcpy.NumPyArrayToRaster(zz,xyOrig,xyCell,xyCell)
grid1.save(os.path.join(arcpy.workspace,outRaster))
strPrj = "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID"\
         "['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],"\
         "UNIT['Degree',0.0174532925199433]]"
arcpy.DefineProjection_management(outRaster,strPrj)
print "Written: %s" % grid1