# # 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