OCEAN file with DMS and CHLO

PYTHON
1

I have a ocean file named ocean.txt that is in the attachment. I used:

to add DMS and CHLO to this file. My new ocean file is in the attachment with name ocean_new.txt.

Then I compare the new ocean file with the ocean file that is existed in the CMAQ v5.5 test case data, and in my new ocean file I have something like this:

float DMS(TSTEP, LAY, ROW, COL) ;
DMS:_FillValue = -999.f ;
DMS:long_name = "DMS " ;
DMS:var_desc = "DMS " ;
DMS:units = "nM " ;
DMS:fill_value = -999.f ;
DMS:grid_mapping = “crs” ;
DMS:missing_value = -999.f ;
DMS:description = "Hansell et al. seawater DMS climatology " ;
float CHLO(TSTEP, LAY, ROW, COL) ;
CHLO:_FillValue = -32767.f ;
CHLO:long_name = "CHLO " ;
CHLO:var_desc = "CHLO " ;
CHLO:units = "mg m^-3 " ;
CHLO:fill_value = -32767.f ;
CHLO:standard_name = “mass_concentration_of_chlorophyll_in_sea_water” ;
CHLO:grid_mapping = “crs” ;
CHLO:missing_value = -32767.f ;
CHLO:reference = “Hu, C., Lee Z., and Franz, B.A. (2012). Chlorophyll-a algorithms for oligotrophic oceans: A novel approach based on three-band reflectance difference, J. Geophys. Res., 117, C01011, doi:10.1029/2011JC007395.” ;
CHLO:display_scale = “log” ;
CHLO:display_min = 0.01f ;
CHLO:display_max = 20.f ;

but in the EPA’s test case ocean file we have:

float DMS(TSTEP, LAY, ROW, COL) ;
DMS:long_name = "DMS " ;
DMS:units = "nM " ;
DMS:var_desc = "DMS " ;
float CHLO(TSTEP, LAY, ROW, COL) ;
CHLO:long_name = "CHLO " ;
CHLO:units = "mg m^-3 " ;
CHLO:var_desc = "Chlorophyll Concentration, OCI Algorithm

I think my new ocean file has some issues since the values of DMS and CHLO in all of the grid cells are equal to 0. In addition we have some additional lines in the variables part of my new ocean file in compare with the EPA’s test case ocean file. I’m not sure why my new ocean file created wrongly.

I will appreciate if you can create a new ocean file with your python code which is available in that mentioned GITHUB website with my attached ocean.txt file.

@barronh

ocean.txt (485.7 KB)
ocean_new.txt (951.2 KB)

2

I ran the Notebook on Google Colab with ocean.txt as an input (after using ncgen) and it worked fine.

image
image1002×628 72.5 KB

Steps

To make this work on Google Colab, you need to upload your ocean file, which I will assume is called “OCEANFile_AIRPACT_04km_fixed.nc”

  1. Follow this link to open on google colab

  2. Drag and drop your ocean file to the file browser.

  3. Add these three cells to the top and run them.

%%writefile requirements.txt
setuptools >= 58.0
importlib_metadata >= 4.6
numpy >= 1.19.5,<2
pandas >= 1.1.5
netCDF4 >= 1.5.8
matplotlib >= 3.3.4
pycno >= 0.2.0
pyproj >= 2.6.1
pseudonetcdf >= 3.2.0
cdo == 1.5.3

!python -m pip install -r requirements.txt

!apt install netcdf-bin
!wget -N https://forum.cmascenter.org/uploads/short-url/unjiMkCmgd6FBTWw3U0SZJ0f6AR.txt
!ncgen -o OCEANFile_AIRPACT_04km_fixed.nc unjiMkCmgd6FBTWw3U0SZJ0f6AR.txt

  1. After running those cells, restart the notebook (Runtime - Restart session)

  2. Update the user specifications

dom = 'AIRPACT_04km'

and

ocnintmpl = '/content/OCEANFile_AIRPACT_04km_fixed.nc'

  1. And, in install prerequisites, set installcdo = True

  2. Run the notebook and download results.

  3. To make a figure like the one I attached, run the code below. The exact result will depend on which file you visualize.

import matplotlib.pyplot as plt
import matplotlib.colors as mc
import pycno  # pip install pycno

tmpf = pnc.pncopen(ocnoutpath, format='ioapi')
cno = pycno.cno(tmpf.getproj(withgrid=True))
fig, axx = plt.subplots(1, 2, figsize=(12, 4))
for ax, key, norm in zip(axx, 'DMS CHLO'.split(), [None, mc.TwoSlopeNorm(4, 0)]):
  Z = tmpf[key]
  qm = ax.pcolormesh(Z[0, 0], norm=norm)
  fig.colorbar(qm, label=f'{Z.long_name} [{Z.units}]'.replace('  ', ''))
cno.drawstates(ax=axx)

p.s., the units on your ocean file are non-ascii characters – not ideal, but not a problem for the notebook.

3 Likes
3

Hi Barron,

Thank you so much for your comment.

I did the same as you mentioned but during the notebook’s run and when I run:

gdf = pnc.pncopen(gdpath, format=‘ioapi’)

proj = gdf.getproj()

gproj = gdf.getproj(withgrid=True)

crs = proj.crs.to_cf()

I get the following error:

CRSError Traceback (most recent call last)
in <cell line: 0>()
1 gdf = pnc.pncopen(gdpath, format=‘ioapi’)
----> 2 proj = gdf.getproj()
3 gproj = gdf.getproj(withgrid=True)
4 crs = proj.crs.to_cf()

5 frames
/usr/local/lib/python3.11/dist-packages/pyproj/_crs.pyx in pyproj._crs._CRS.init()

CRSError: Invalid projection: +proj=lcc +lat_1=np.float64(30.0) +lat_2=np.float64(60.0) +lat_0=np.float64(49.0) +lon_0=np.float64(-121.0) +y_0=np.float64(942000.0) +x_0=np.float64(342000.0) +a=6370000.0 +b=6370000.0 +no_defs +type=crs: (Internal Proj Error: proj_create: Error 1027 (Invalid value for an argument): lcc: Invalid value for lat_1 and lat_2: |lat_1 + lat_2| should be > 0)

I will appreciate it if you have time to check it.
@barronh

4

Thanks for the extra detail. That’s a problem with the numpy version 2 that I am addressing in a new version of PseudoNetCDF for now, downgrade numpy. I edited the instructions above to ensure that it works.

1 Like
5

Thanks Barron.

It worked.