From:	"Jason Roberts" <>
To:	"'Tim Sippel'" <>, <>
Subject:	RE: [mget-help] extract netcdf 3D
Date:	Mon, 13 Jul 2009 20:35:30 -0400

Tim,   The short story is, no, we haven't made any progress on that yet, but that it might be available sometime in September.   The longer story is that we are planning to develop a set of modules that will simplify conversion and sampling of gridded data from a variety of formats, as part of a project that starts later this summer. That project's specific scenario involves sampling 3D netCDFs (the ROMS-CoSINE Pacific ocean model) so it is likely that we'll develop a 3D netCDF tool for it.   In the mean time, if you are interested in writing your own Python script, you could use the attached one as an example. This script builds monthly and cumulative climatologies from ROMS-CoSINE netCDFs that I received from Fei Chai, one of the authors of that model. Each netCDF contained one year of data. A given variable was 3D, with 12 time slices, one for each month of that year. The script takes four arguments: a "glob" _expression_ specifying the netCDFs to process, a variable name from the netCDF to process, and output monthly and cumulative climatology rasters to create.   The script uses Jeff Whittaker's excellent netcdf4-python library to read netCDFs, and numpy to manipulate them. Copies of both of these are assimilated into MGET, so MGET can access them without users performing an additional install.   Hope that helps,   Jason   From: Tim Sippel [mailto:] Sent: Monday, July 13, 2009 6:35 PM To: Subject: [mget-help] extract netcdf 3D   Hi Jason- I just read an email thread from a few months back where you mention developing an MGET tool to extract netcdf 3D data.  Has this been implemented yet?  Thanks, Tim

import glob
import os
import sys

from GeoEco.ArcGIS import GeoprocessorManager
from GeoEco.DataManagement.ArcGISRasters import ArcGISRaster
from GeoEco.Dependencies import PythonAggregatedModuleDependency
import GeoEco.AssimilatedModules.netCDF3

# Parse command line arguments.

globExpr = sys.argv[1]
varName = sys.argv[2]
monthlyRaster = sys.argv[3]
cumulativeRaster = sys.argv[4]

# Initialize numpy from MGET and import it.

d = PythonAggregatedModuleDependency('numpy')
d.Initialize()
import numpy

# Read the netcdfs into numpy arrays.

monthlyArrays = None
cumulativeArray = None
totalMonths = 0.0

netCDFs = glob.glob(globExpr)
for netCDF in netCDFs:
    print('Reading %s.' % netCDF)
    dataset = GeoEco.AssimilatedModules.netCDF3.Dataset(netCDF, 'r')
    
    if monthlyArrays is None:
        monthlyArrays = []
        for i in range(12):
            
monthlyArrays.append(numpy.zeros(dataset.variables[varName].shape[2:], 
'float32'))
        cumulativeArray = numpy.zeros(dataset.variables[varName].shape[2:], 
'float32')
    
    for i in range(12):
        data = dataset.variables[varName][i,0,:,:]
        monthlyArrays[i] += data
        cumulativeArray += data
        totalMonths += 1.0

# Flip the arrays and divide by the total number of months.

for i in range(12):
    monthlyArrays[i] = numpy.flipud(monthlyArrays[i]) / 
(numpy.zeros(monthlyArrays[i].shape, 'float32') + totalMonths)

cumulativeArray = numpy.flipud(cumulativeArray) / 
(numpy.zeros(cumulativeArray.shape, 'float32') + totalMonths)

# Write the rasters.

GeoprocessorManager.InitializeGeoprocessor()

print('Writing %s.' % cumulativeRaster)
ArcGISRaster.FromNumpyArray(cumulativeArray, cumulativeRaster, 
-8907324.838256, -5595576.075406, 13880, nodataValue=0, 
coordinateSystem="PROJCS['World_Mercator',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-180.0],PARAMETER['Standard_Parallel_1',0.0],UNIT['Meter',1.0]]
 # 
GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
 buildPyramids=True, overwriteExisting=True)

for i in range(12):
    print('Writing %s%02i.' % (monthlyRaster, i+1))
    ArcGISRaster.FromNumpyArray(monthlyArrays[i], '%s%02i' % (monthlyRaster, 
i+1), -8907324.838256, -5595576.075406, 13880, nodataValue=0, 
coordinateSystem="PROJCS['World_Mercator',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-180.0],PARAMETER['Standard_Parallel_1',0.0],UNIT['Meter',1.0]]
 # 
GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]",
 buildPyramids=True, overwriteExisting=True)

< date >
< thread >

Archives powered by MHonArc.