Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

UPDATE: 7/26/2013

The issue I'm encountering is a known bug with the netCDF-Java API described in ticket TDS-391:

A user would like to convert radar data to netCDF format. Unfortunately, the data files have structure data, which fails badly for writing netcdf-3 (illegal data type, as netcdf-3 does not support structures) or netcdf-4 (unimplemented) files.

The issue is targeted for a fix in API ver. 4.4, which will hopefully will be out this summer, as per an email from Unidata netCDF Java Support received on July 8:

We are going to implement a solution for this in version 4.4, but it won't be out before our training workshop later this month. We may limit the implementation to NetCDF-4, and offer to skip trying to write structures in NetCDF-3.


I have many level-III NEXRAD precipitation netCDF files which I created from the raw files using the netCDF-Java API (ver. 4.3.17). I am now trying to compute monthly average rainfall rates and returning a new netCDF file which I will then need to convert to ESRI GRID format. I am new to working with netCDF files, but have read up a bit on creating dimensions and variables using the scipy.io.netcdf module.

I am using the first netCDF file in the directory as a sort of template file from which I'm grabbing the x and y locations and projection information. I have 3 dimensions:

In [4]: ncfile.dimensions
Out[4]: {'textStringSize1': 1, 'x': 131, 'y': 131}

I also have 4 variables:

In [3]: ncfile.variables
Out[3]: 
{'FlatEarth': <scipy.io.netcdf.netcdf_variable at 0x67527d0>,
 'PrecipArray_0': <scipy.io.netcdf.netcdf_variable at 0x6752690>,
 'x': <scipy.io.netcdf.netcdf_variable at 0x67526d0>,
 'y': <scipy.io.netcdf.netcdf_variable at 0x6752430>}

I have successfully created my x and y dimensions and variables and populated the precip_array_0 variable with my monthly average values. The last piece of the puzzle, I think, is to copy the projection info which I believe is what the FlatEarth variable holds. I'm not sure how to look at what's in the variable:

In [5]: ncfile.variables['FlatEarth'][:]
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
C:\Users\jbellino\Desktop\<ipython-input-5-aec6027d7a77> in <module>()
----> 1 ncfile.variables['FlatEarth'][:]

C:\Python27_epd32\lib\site-packages\scipy\io\netcdf.pyc in __getitem__(self, index)
    837 
    838     def __getitem__(self, index):
--> 839         return self.data[index]
    840 
    841     def __setitem__(self, index, data):

IndexError: 0-d arrays can't be indexed.


In [6]: ncfile.variables['FlatEarth']
Out[6]: <scipy.io.netcdf.netcdf_variable at 0x67527d0>


In [10]: ncfile.variables['FlatEarth'][()]
Out[10]: ' '

Calling help on the object returns help for the entire variables object, not just the FlatEarth variable:

In [11]: ncfile.variables['FlatEarth']??
Type:       dict
Base Class: <type 'dict'>
String Form:{'y': <scipy.io.netcdf.netcdf_variable object at 0x06752430>, 'x': <scipy.io.netcdf.netcdf_variab <...> e object at 0x067527D0>, 'PrecipArray_0': <scipy.io.netcdf.netcdf_variable object at 0x06752690>}
Namespace:  Interactive
Length:     4
Docstring:
dict() -> new empty dictionary
dict(mapping) -> new dictionary initialized from a mapping object's
    (key, value) pairs
dict(iterable) -> new dictionary initialized as if via:
    d = {}
    for k, v in iterable:
        d[k] = v
dict(**kwargs) -> new dictionary initialized with the name=value pairs
    in the keyword argument list.  For example:  dict(one=1, two=2)

I'm assuming that the FlatEarth variable uses the textStringSize1 dimension, but I'm not sure how to create this variable in the new netCDF file or populate it once it's been created. Any help is much appreciated!

share|improve this question

2 Answers 2

You could avoid writing the NetCDF file altogether and just use arcpy.NumPyArrayToRaster to bring your numpy array of mean precip directly into ArcMap. The code could look someting like this:

import arcpy
xyOrig = arcpy.Point(float(lonmin),float(latmin))
arcpy.workspace  = "C:\\workspace"
arcpy.env.overwriteOutput = True
rasterName = "traster"
outRaster = os.path.normpath(os.path.join(arcpy.workspace,rasterName))
grid1=arcpy.NumPyArrayToRaster(z,xyOrig,dx,dy)
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
arcpy.AddMessage("Written: %s" % grid1)

The full example from which this snippet was obtained is here:
https://github.com/rsignell-usgs/dap2arc/blob/master/dap2raster.py

BTW, when I need to write NetCDF, I use the NetCDF4-Python library as it can read and write NetCDF4 files and read from OPeNDAP in addition to reading and writing NetCDF3 files. It's nice to use a common interface for all. And you can get a 32-bit windows version of NetCDF4-Python compatible with ArcGIS 10.0+ at Gohlke's Python Binaries for Windows site

Good luck!

share|improve this answer
    
Hi Rich, thanks for the code. The only problem with going straight to ESRI grid is that my x,y locations the point of origin is unknown (e.g., 0,0; 4,4; 20,92, etc.). I have no way of placing these in a GCS without the projection info. I'll try the netcdf4 package from Gohlke's site to see if there are additional bits of information that I'm not getting with the scipy module. –  Jason Jun 17 '13 at 13:35
    
I don't quite understand. Do you mean you don't know what the coordinate system of your x,y points are? –  Rich Signell Jun 18 '13 at 22:03
    
Yes, I can't access the coordinate system information embedded within the netCDF (or it's not getting translated from the raw nexrad files). I don't know where the origin lies in space. I'm starting to go down a rabbit hole with this dataset - I'm trying to compare computed effective precipitation using the nexrad and prism datasets, but time on the project is running out. Not a huge deal if I can't get this worked out. –  Jason Jun 19 '13 at 12:54
    
It looks like the NetCDF Java library I'm using to decode the raw nexrad files is not working correctly. If/when I get things squared away with the decoding process I will revisit this page to update. –  Jason Jun 19 '13 at 16:12

From my limited experience working with NetCDF files, coordinate systems can still be problematic. They are supported in extensions like the CF conventeions, and this section of the CF Metadata standard on coordinate reference systems is a good place to start. I believe that the reference to flat_earth means that your netCDF files are intentionally ignoring the sphere, and the details for how the variable is specified is included within the NEXRAD documentation. As a fallback, perhaps you can copy this information directly?

I'm not sure of your entire pipeline, but if you are ending up with ESRI GRID files, you may want to consider injecting the projection information into the data at that point, which is probably easier than getting it correctly embedded in the NetCDF file itself.

share|improve this answer
    
Thanks for the links. –  Jason Jun 12 '13 at 20:24
    
I'm not entirely sure how to get the projection information from the netCDF files to define the projection for the resulting rasters. The FlatEarth variable value appears to be " " - this is the value returned when I call ncfile.variables['FlatEarth'][()]. The NOAA Weather and Climate Toolkit appears to do a bit of magic when you convert the raw NEXRAD files using that GUI - it converts 0-255 values to hourly rates of rainfall (for the digital precip product I have). Maybe it also applies the projection information as well? –  Jason Jun 13 '13 at 12:43

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.