Not logged inOcean Color Forum

The forum is locked.

The Ocean Color Forum has transitioned over to the Earthdata Forum (https://forum.earthdata.nasa.gov/). The information existing below will be retained for historical reference. Please sign into the Earthdata Forum for active user support.

Up Topic Products and Algorithms / Satellite Data Products & Algorithms / adding a product to L2 in Python (locked)
- By peland Date 2016-03-30 15:52
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
- By gnwiii Date 2016-03-30 17:06
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'),
- By peland Date 2016-03-31 16:08
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
- By melliott Date 2016-03-31 17:43
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
- By gnwiii Date 2016-03-31 17:57
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)
Up Topic Products and Algorithms / Satellite Data Products & Algorithms / adding a product to L2 in Python (locked)

Powered by mwForum 2.29.7 © 1999-2015 Markus Wichitill