barronh / pseudonetcdf

PseudoNetCDF like NetCDF except for many scientific format backends
GNU Lesser General Public License v3.0
76 stars 35 forks source link

Converting WRF NetCDF to NetCDF-CF format #109

Closed nishadhka closed 4 years ago

nishadhka commented 4 years ago

Many thanks for library.

I am trying to convert WRF output in NetCDF into NetCDF-CF format to get projected (x.y) array for certain variables.

The sample WRF output is attached here with https://drive.google.com/file/d/1rg9IG_dJG_Iz5kKOFBXSf6K3fSx7p4n1/view?usp=sharing

1.using Python import

import PseudoNetCDF as pnc

inncfile='wrfout_d03_2020-07-23_00:00:00'
outputfile='wrfout_cf.nc'
infile = pnc.pncopen(inncfile, addcf=True)
infile.save(outputfile)

Which ends up in error of

>>> infile.save(outputfile)
Adding dimensions
Adding globals
Adding variables
Defining Times
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/core/_files.py", line 2329, in save
    return pncwrite(self, *args, **kwds)
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/_getwriter.py", line 42, in pncwrite
    return pncgen(*args, **kwds)
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncgen.py", line 245, in pncgen
    format=format)
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncgen.py", line 52, in convert
    self.addVariables(pfile, nfile)
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncgen.py", line 165, in addVariables
    self.addVariable(pfile, nfile, k, data=self.datafirst)
  File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncgen.py", line 134, in addVariable
    k, typecode, pvar.dimensions, **create_variable_kwds)
  File "netCDF4/_netCDF4.pyx", line 2808, in netCDF4._netCDF4.Dataset.createVariable
  File "netCDF4/_netCDF4.pyx", line 3804, in netCDF4._netCDF4.Variable.__init__
TypeError: illegal primitive data type, must be one of dict_keys(['S1', 'i1', 'u1', 'i2', 'u2', 'i4', 'u4', 'i8', 'u8', 'f4', 'f8']), got |S0
  1. used command line tool as follows
    
    $$$ pncgen.py -f uamiv wrfout_d03_2020-07-23_00:00:00 test.nc
    Traceback (most recent call last):
    File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncparse.py", line 895, in getfiles
    f = pncopen(ipath, format=file_format, **format_options)
    File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/_getreader.py", line 153, in pncopen
    outfile = reader(*args, **kwds)
    File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 131, in __init__
    self.__readheader()
    File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 263, in __readheader
    GDTYPE = self.GDTYP = {0: 1, 1: 5, 2: 2, 3: 6}[cproj]
    KeyError: 1413808128

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncparse.py", line 901, in getfiles allreaders[file_format](ipath, **format_options) File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 131, in init self.readheader() File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/camxfiles/uamiv/Memmap.py", line 263, in readheader GDTYPE = self.GDTYP = {0: 1, 1: 5, 2: 2, 3: 6}[cproj] KeyError: 1413808128

During handling of the above exception, another exception occurred:

Traceback (most recent call last): File "/home/Nishadh/miniconda3/bin/pncgen.py", line 3, in main() File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncgen.py", line 269, in main ifiles, options = pncparse(has_ofile=True, interactive=False) File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncparse.py", line 611, in pncparse ifiles.extend(pncprep(subarg)[0]) File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncparse.py", line 839, in pncprep fs = getfiles(ipaths, args) File "/home/Nishadh/miniconda3/lib/python3.7/site-packages/PseudoNetCDF/pncparse.py", line 911, in getfiles raise oute # from e OSError: Unable to open path with uamiv(path, **{}) path="wrfout_d03_2020-07-23_00:00:00" error="1413808128"


Any pointers to address this issue will be a great help, thanks. 
barronh commented 4 years ago

There are three separate issues, and I will explain each. None of these problems are show stoppers, depending on what you want to do. I'll show a working example and you can decide if you want to pursue this in PseudoNetCDF or wrf-python (See link below).

First, The support for WRF meta-data is somewhat limited. In my experience, the time coordinate as text is sensitive to versions of numpy and the netcdf reader. So, reading is easy but writing it is harder. Further, WRF just hasn't been extensively used with PseudoNetCDF, so several things are harder than they need to be. Depending on what you want to do, this may be a trivial barrier. It is also an opportunity for me to make some updates.

Second, the uamiv example you're using is trying to read a wrf netcdf file as if it was a CAMx unformatted fortran file. That won't work. The -f option from the command line tells PseudoNetCDF what kind of file it is opening. The uamiv is the Urban Air Shed Model version IV, which CAMx adopted with adaptations. So, PseudoNetCDF is reading the data as if it was uamiv and it fails because it is not.

Third, the addcf option really only works for uamiv and ioapi formats. I fails silently, which I should update.

Now, depending on what you want to do, this may not be a problem. Below is an example that plots T2 in degrees celcius.

import PseudoNetCDF as pnc

wrfpath = 'wrfout_d03_2020-07-23_00_00_00'

rawf = pnc.pncopen(wrfpath, format='netcdf')
cff = rawf.sliceDimensions(bottom_top=0)

# Basically addcf asa a function.
pnc.conventions.ioapi.add_cf_from_wrfioapi(cff)

# Adds meta variables: time, longitude, latitude, x, y,
# and lambert_conformal_conic
# x, y are projected units, but no not include the false_easting/northing
# lambert_conformal_conic is CF projection definition including false_easting/northing
# in the future, x/y west_east, south_north should be added to match the lambert_conformal
# the time variable is not correct, so be careful.
# for now, I add these variables manually
x = cff.variables['x']
y = cff.variables['y']
cff.copyVariable(x - x.min(), key='west_east')
cff.copyVariable(y - y.min(), key='south_north')

# Using eval to convert temperature to celcius. This is a trivial example, but useful
# for computing WRF derived variables.
cff.eval('T2C = T2 - 273.15; T2C.units = "C"', inplace=True)

# Now make a plot using the coordinates added above
ax = cff.plot('T2C', plottype='west_east-south_north')

# Adding a map is usually part of the plot function, but needs to 
# be added manually via the getMap function.
bmap = cff.getMap(
    maptype='basemap',
    llcrnrx=x.min(), urcrnrx=x.max(),
    llcrnry=y.min(), urcrnry=y.max(),
)
bmap.drawcoastlines(ax=ax)
bmap.drawcountries(ax=ax)
# Selected some parallels and meridians based on yoru data
bmap.drawparallels([5, 10], labels=['5N', '10N'])
bmap.drawmeridians([71, 72], labels=['71E', '72E']);

# Save as a figure.
ax.figure.savefig('WRF_T2C.png')

del cff.variables['Times']
cff.save('test.nc')

Like I said, I don't know what you want to do so this is just a toy example that may or may not help. See the link below for an alternative system.

Link to WRF example with xarray: https://scott-hosking.github.io/notebooks/wrf_with_xarray/

nishadhka commented 4 years ago

Many thanks for the reply. The sample code works perfectly.

Sorry missed to explain the need for conversion. The aim is to reproject WRF output in NetCDF form into Geographic coordinate system. The resultant output, for example, a single time T2 which can be represented as 2D array with its geographic information represented as [start_longitude, start_latitude, no_of_xaxis_grid, no_of_yaxis_grid, grid_resolution].

There are issues related to this conversion as described in this post https://fabienmaussion.info/2018/01/06/wrf-projection/. So far tried with http://www.meteo.unican.es/wiki/cordexwrf/SoftwareTools/WrfncXnj, but this is Python2 version and wgrib2 https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/ has facility but the input data I guess needs to be in grib form.

Had a trail with wrf-python and needed to explore on its interpolation methods.

Following the given code in your earlier message for projecting the WRF output in LCC. The steps

x = cff.variables['x']
y = cff.variables['y']
cff.copyVariable(x - x.min(), key='west_east')
cff.copyVariable(y - y.min(), key='south_north')

Does generate the x and y in lcc form

import xarray as xr

db=xr.open_dataset('test.nc')
x_lat=db['x'].values
x_lon=db['y'].values
x_lat, x_lon
array([-1048001.54786551, -1044001.54786551, -1040001.54786551,
        -1036001.54786551, -1032001.54786551, -1028001.54786551,
        -1024001.54786551,....]
array([975999.42921298, -971999.42921298, -967999.42921298,
        -963999.42921298, -959999.42921298, -955999.42921298,
        -951999.4292129....]

Is there any step related to convert the lcc into geographic coordinate(lat,lon) system. My current plan is use custom converter as described in the link https://stackoverflow.com/a/9898274/2501953.

barronh commented 4 years ago

Several things to note. First, WRF gave you latitude and longitude as XLAT and XLONG. Second, those variables are already renamed as latitude/longitude by by add_cf_from_wrfioapi. The add_cf_from_wrfioapi just renames them. So, WRF gives you the lat/lon coordinates.

So, you have the coordinates if you need them. No problem there. PseudoNetCDF provides convenience functions to convert between i,j and x,y and lon,lat (cff.ll2ij, ij2ll, ll2xy, xy2ll, etc). So, getting the coordinates is no problem.

If you actually want to "reproject", you can combine the PsedoNetCDF code with cdo. The two code blocks below do just that. The first makes a cdo-compliant file (test.nc). The second uses cdo remapbil to bilinearly interpolate to a lon/lat grid (wrflonlat.nc).

These python commands below produce a file with just U, V, and T2 that can be used with cdo to remap.

import numpy as np
import PseudoNetCDF as pnc

wrfpath = 'wrfout_d03_2020-07-23_00_00_00'

rawf = pnc.pncopen(wrfpath, format='netcdf')
cff = rawf.sliceDimensions(bottom_top=0)

# Basically addcf asa a function.
pnc.conventions.ioapi.add_cf_from_wrfioapi(cff)

# Adds meta variables: time, longitude, latitude, x, y,
# and lambert_conformal_conic
# x, y are projected units, but no not include the false_easting/northing
# lambert_conformal_conic is CF projection definition including false_easting/northing
# in the future, x/y west_east, south_north should be added to match the lambert_conformal
# the time variable is not correct, so be careful.
# for now, I add these variables manually
x = cff.variables['x']
y = cff.variables['y']
cff.copyVariable(x - x.min() + cff.DX / 2, key='west_east')
cff.copyVariable(y - y.min() + cff.DY / 2, key='south_north')
xs = cff.copyVariable(x, key='west_east_stag', dimensions=('west_east_stag',), withdata=False)
xs[:] = np.arange(len(cff.dimensions['west_east_stag'])) * cff.DX
ys = cff.copyVariable(y, key='south_north_stag', dimensions=('south_north_stag',), withdata=False)
ys[:] = np.arange(len(cff.dimensions['south_north_stag'])) * cff.DY

del cff.variables['Times']
# delattr(cff.variables['west_east'], 'coordinates')
# delattr(cff.variables['south_north'], 'coordinates')
# delattr(cff.variables['T2'], 'coordinates')
cff.subsetVariables(['lambert_conformal_conic', 'west_east', 'south_north', 'west_east_stag', 'south_north_stag', 'U', 'V', 'T2']).save('test.nc')

The cdo commands below (adapted from https://code.mpimet.mpg.de/boards/2/topics/96) would use bilinear interpolation to put on a lonlat grid.

cat > mygrid << EOF
gridtype = lonlat
xsize    = 360
ysize    = 180
xfirst   = -179.5
xinc     = 1
yfirst   = -89.5
yinc     = 1
EOF

cdo remapbil,mygrid test.nc wrflonlat.nc

You can edit the my grid file to produce a more appropriate grid. This grid is good for testing because it is fast and easy to confirm on a map. Note also, you could choose to use other remap functions (e.g., remapycon, remapnn, etc).

As a final note, some projections require newish versions of cdo. For LCC, I am not sure it matters as much. In this example, I was using cdo v1.9.6.

nishadhka commented 4 years ago

Had a test run with Pseudonetcdf and CDO. CDO is new to me and thanks for introducing. Facing some issue with opening the netcdf file generated from Pseudonetcdf step mentioned in earlier step. The cdo info command ends up in this error. This error I also had earlier. as it is randomly popping up and even observed for reading the output file from xarray lib. I am looking more into this error, will raise a new issue related to this. Used Anaconda provided cdo, https://anaconda.org/conda-forge/cdo, used version 1.9.6 instead of latest as well but got this error.

$$ cdo -info test.nc
HDF5-DIAG: Error detected in HDF5 (1.10.4) thread 139821330716480:
  #000: H5F.c line 509 in H5Fopen(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #001: H5Fint.c line 1400 in H5F__open(): unable to open file
    major: File accessibilty
    minor: Unable to open file
  #002: H5Fint.c line 1615 in H5F_open(): unable to lock the file
    major: File accessibilty
    minor: Unable to open file
  #003: H5FD.c line 1640 in H5FD_lock(): driver lock request failed
    major: Virtual File Layer
    minor: Can't update object
  #004: H5FDsec2.c line 941 in H5FD_sec2_lock(): unable to lock file, errno = 11, error message = 'Resource temporarily unavailable'
    major: File accessibilty
    minor: Bad file ID accessed

cdo info: Open failed on >test.nc<
Unknown Error

Closing this issue as of now.