Hi All,

I'm trying to add a modified chlorophyll product to a MODIS L2 file. Here's my code:

import numpy as np,shutil
from netCDF4 import Dataset
shutil.copy(inFile, outFile)
nc=Dataset(outFile, 'a')
gd=nc.groups['geophysical_data']
v0=gd.variables['chlor_a']
v=gd.createVariable('chl2','f',('number_of_lines','pixels_per_line'),
fill_value=v0._FillValue)
v.long_name='Corrected '+v0.long_name
v.units=v0.units
v.standard_name=v0.standard_name
v.valid_min=v0.valid_min
v.valid_max=v0.valid_max
v.reference=v0.reference
v[:]=chl2
nc.close()

The resulting file looks good in Python, but when I try to open it in SeaDAS7.3, I get

An exception occurred:
Type: java.util.concurrent.ExecutionException
Message: java.lang.ArrayIndexOutOfBoundsException: -1354

So it looks like something in the netcdf still thinks it has the old number of variables. Can you help?

Cheers, Peter

You used:

v=gd.createVariable('chl2','f',('number_of_lines','pixels_per_line'),


I've done something similar, but using R, which gives me:


float geophysical_data/chlor_a[pixels_per_line,number_of_lines]   (Chunking: [1067,10])  (Compression: level 4)"


Unless Python reverses the order of array indices, you could try:


v=gd.createVariable('chl2','f',('pixels_per_line','number_of_lines'),


Hi,

Thanks for the tip! Would you mind showing me your R code? I'm using R to do stats then switching to Python for the netCDF manipulation, so it would help to be able to do the netCDF stuff in R.

Cheers, Peter

Peter,

What is chl2 supposed to hold? It isn't defined in the code you posted. When I tried to run your code, it produced an error because the chl2 variable used in the next to last line isn't defined. An output file was produced, but it didn't have the chl2 you were looking for.

Just to make the code run, I commented out the offending line. That produced an output file that, when opened in SeaDAS, did have chl2, although it held no data. I think if you set up that chl2 variable you'll get what you're looking for.

HTH

-Matt

I was wrong about needing to reverse the indices.  Here is an example to create a band with (chlor_a - chl_ocx):


inFile='../V2015127150519.L2_NPP_OC.nc'
outFile='V2015127150519.L2_NPP_OC.nc'
# http://oceancolor.gsfc.nasa.gov/forum/oceancolor/topic_show.pl?pid=26505
import numpy as np 
import shutil
from netCDF4 import Dataset
shutil.copy(inFile, outFile)
nc=Dataset(outFile, 'a')
gd=nc.groups['geophysical_data']
chlor_a=gd.variables['chlor_a']
print(chlor_a)
chl_ocx=gd.variables['chl_ocx']
v=gd.createVariable('chl_dif','f',('number_of_lines','pixels_per_line'),
   fill_value=chlor_a._FillValue)
v.long_name=chlor_a.long_name + ' - ' + chl_ocx.long_name
v.units=chlor_a.units
v.standard_name='chl_dif'
v.valid_min=-0.1
v.valid_max=0.1
v.reference=chlor_a.reference
v[:]=chlor_a[:]-chl_ocx[:]
print(v)
nc.close()


I did discover a while back that I needed a current Unidata netcdf4 Python version.  I have:




>>> import netCDF4
>>> print(netCDF4.__version__)

1.2.1


Here is the R code:

require(ncdf4)
l2dir="<level2 directory>"
l2ncfile='<filename>'
# Add a new bands for the difference and ratio between `chlor_a` and `chl_ocx`:
nc = nc_open(file.path(l2dir,l2ncfile), write=TRUE)
chlor_a = ncvar_get(nc, 'geophysical_data/chlor_a')
chl_ocx = ncvar_get(nc, 'geophysical_data/chl_ocx')
chl_dif = chlor_a - chl_ocx
chl_rat = chl_ocx/chlor_a

pixels_per_line = nc$dim[["pixels_per_line"]]
number_of_lines = nc$dim[["number_of_lines"]]

var_chl_dif = ncvar_def('geophysical_data/chl_dif', 'mg m^-3', list(pixels_per_line, number_of_lines), NA, "Chlorophyll Concentration Difference, OCI-OCX")
var_chl_rat = ncvar_def('geophysical_data/chl_rat', 'mg m^-3', list(pixels_per_line, number_of_lines), NA, "Chlorophyll Concentration Ratio, OCX/OCI")
nc = ncvar_add(nc, var_chl_dif)
nc = ncvar_add(nc, var_chl_rat)
ncvar_put(nc, var_chl_dif, chl_dif)
ncatt_put(nc, 'geophysical_data/chl_dif', "valid_min", -0.1, prec="float")
ncatt_put(nc, 'geophysical_data/chl_dif', "valid_max", 0.1, prec="float")
ncvar_put(nc, var_chl_rat, chl_rat)
ncatt_put(nc, 'geophysical_data/chl_rat', "valid_min", 0)
ncatt_put(nc, 'geophysical_data/chl_rat', "valid_max", 10)
nc_close(nc)