Add Support for Reading/Writing Intermediate WRF Files

[bladwig1]

Had a user request the ability to read/write the intermediate files. Although most users use WPS, should probably support the WRFSI and MM5 as well. See

http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Basics/IM_files/

[lesserwhirls]

HI @bladwig1 - are you still looking into writing intermediate files?

[bladwig1]

@lesserwhirls - This is still planned, but it won't be in the soon-to-be-released 1.1. I do have a really primitive reader for WPS format (I don't know anyone that uses WRFSI or MM5 format). If you need something to hack up, I can send you the reader and you can reverse it to do the writer.

[lesserwhirls]

What I'd like to do is enable the use of the the NetcdfSubset Service (NCSS) on the THREDDS Data Server (TDS) in the WRF initialization process. While the TDS has no problem serving out GRIB data, NCSS only returns netCDF files, which is somewhat of a no-go for WRF initialization at the moment.

We've danced around the idea of extending ungrib to handle netCDF, but we at Unidata are no experts in WRF and would need a bit of help there. As a way around that, we've looked at the possibility of writing a netCDF -> intermediate file converter, driven by WRF vtables. I wanted to check in on this issue to see if you had already started on the path of making an intermediate format writer yet. If not, I will tackle that in the next month or so, and would be happy to contribute that back to this project (if interested).

[bladwig1]

I sent an email with the WPS Python reader and also the Fortran code that NCL uses. Hopefully it will save you some time. Feel free to contribute back what you have done, since I probably won't get to this task for a while.

[lesserwhirls]

Awesome @bladwig1 - thank you!

[DonteRW]

@bladwig1 @lesserwhirls

hie

any progress on this issue?

How can i get this "WPS Python reader and also the Fortran code that NCL uses"?

[bladwig1]

The documentation for the intermediate file format is below. I think NCL only supports WPS format. I'm not sure if anyone still uses MM5 or WRFSI format.

http://www2.mmm.ucar.edu/wrf/OnLineTutorial/Basics/IM_files/

The fortran code from NCL can be found at these links:

https://github.com/NCAR/ncl/blob/master/ni/src/lib/nfpfort/plotfmt_open.f90 https://github.com/NCAR/ncl/blob/master/ni/src/lib/nfpfort/plotfmt_rddata.f https://github.com/NCAR/ncl/blob/master/ni/src/lib/nfpfort/plotfmt_rdhead.f90 https://github.com/NCAR/ncl/blob/master/ni/src/lib/nfpfort/wrf_write_wps.f

The clunky Python WRF intermediate reader is below. You'll probably want to rewrite it.

from __future__ import print_function
import struct

import numpy as np

# Number of bytes used at the start and end of a fortran record to 
# indicate the record size
SIZE_BYTES = 4

class WPSData(object):
    def __init__(self, ifv=None, hdate=None, xfcst=None, map_source=None, 
                 field=None, units=None, desc=None, xlvl=None, nx=None, 
                 ny=None, iproj=None, startloc=None, startlat=None, 
                 startlon=None, deltalat=None, deltalon=None, nlats=None, 
                 dx=None, dy=None, xlonc=None, truelat1=None, truelat2=None, 
                 earth_radius=None, is_wind_earth_rel=None, slab=None):

        self.ifv = ifv
        self.hdate = hdate
        self.xfcst = xfcst
        self.map_source = map_source
        self.field = field
        self.units = units
        self.desc = desc
        self.xlvl = xlvl
        self.nx = nx
        self.ny = ny
        self.iproj = iproj
        self.startloc = startloc
        self.startlat = startlat
        self.startlon = startlon
        self.deltalat = deltalat
        self.deltalon = deltalon
        self.nlats = nlats
        self.dx = dx
        self.dy = dy
        self.xlonc = xlonc
        self.truelat1 = truelat1
        self.truelat2 = truelat2
        self.earth_radius = earth_radius
        self.is_wind_earth_rel = is_wind_earth_rel
        self.slab = slab

def _parse_record1(data):
    result = {}
    result["ifv"] = struct.unpack(">i", data)

    return result

def _parse_record2(data):
    result = {}
    parsed = struct.unpack(">24sf32s9s25s46sfiii", data)
    result["hdate"] = parsed[0]
    result["xfcst"] = parsed[1]
    result["map_source"] = parsed[2]
    result["field"] = parsed[3]
    result["units"] = parsed[4]
    result["desc"] = parsed[5]
    result["xlvl"] = parsed[6]
    result["nx"] = parsed[7]
    result["ny"] = parsed[8]
    result["iproj"] = parsed[9]

    return result

def _parse_record3(data, iproj):
    result = {}
    if iproj == 0:
        fmt = ">8sfffff"
        parsed = struct.unpack(fmt, data)
        result["startloc"] = parsed[0]
        result["startlat"] = parsed[1]
        result["startlon"] = parsed[2]
        result["deltalat"] = parsed[3]
        result["deltalon"] = parsed[4]
        result["earth_radius"] = parsed[5]
    elif iproj == 1:
        fmt = ">8sffffff"
        parsed = struct.unpack(fmt, data)
        result["startloc"] = parsed[0]
        result["startlat"] = parsed[1]
        result["startlon"] = parsed[2]
        result["dx"] = parsed[3]
        result["dy"] = parsed[4]
        result["truelat1"] = parsed[5]
        result["earth_radius"] = parsed[6]
    elif iproj == 3:
        fmt = ">8sffffffff"
        parsed = struct.unpack(fmt, data)
        result["startloc"] = parsed[0]
        result["startlat"] = parsed[1]
        result["startlon"] = parsed[2]
        result["dx"] = parsed[3]
        result["dy"] = parsed[4]
        result["xlonc"] = parsed[5]
        result["truelat1"] = parsed[6]
        result["truelat2"] = parsed[7]
        result["earth_radius"] = parsed[8]
    elif iproj == 4:
        fmt = ">8sfffff"
        parsed = struct.unpack(fmt, data)
        result["startloc"] = parsed[0]
        result["startlat"] = parsed[1]
        result["startlon"] = parsed[2]
        result["nlats"] = parsed[3]
        result["deltalon"] = parsed[4]
        result["earth_radius"] = parsed[5]
    elif iproj == 5:
        fmt = ">8sfffffff"
        parsed = struct.unpack(fmt, data)
        result["startloc"] = parsed[0]
        result["startlat"] = parsed[1]
        result["startlon"] = parsed[2]
        result["dx"] = parsed[3]
        result["dy"] = parsed[4]
        result["xlonc"] = parsed[5]
        result["truelat1"] = parsed[6]
        result["earth_radius"] = parsed[7]

    return result

def _parse_record4(data):
    result = {}
    result["is_wind_earth_rel"] = struct.unpack(">i", data)

    return result

def _parse_record5(data, nx, ny):
    result = {}

    size = nx * ny
    fmt = ">{}f".format(size)
    parsed = struct.unpack(fmt, data)

    arr = np.array(parsed, dtype=np.float32)
    result["slab"] = arr.reshape((nx, ny), order="F")

    return result

_PARSE_MAP = {0 : _parse_record1,
             1 : _parse_record2,
             2 : _parse_record3,
             3 : _parse_record4,
             4 : _parse_record5}

def parse_record(record_idx, data, iproj=None, nx=None, ny=None):

    if record_idx == 0 or record_idx == 1 or record_idx == 3:
        return _PARSE_MAP[record_idx](data)
    elif record_idx == 2:
        return _PARSE_MAP[record_idx](data, iproj)
    elif record_idx == 4:
        return _PARSE_MAP[record_idx](data, nx, ny)

def read_wps(wps_file, field=None):
    with open(wps_file, "rb") as f:
        wps_params = {}
        record_list = []
        raw_data = f.read()

        record_size_idx = 0
        end_of_file_idx = len(raw_data) - 1

        while True:
            iproj = None
            nx = None
            ny = None
            keep_record = True
            for record_idx in range(5):
                # Each record has the size (in SIZE_BYTES bytes) at the 
                # start and end of each record.  This might be compiler 
                # dependent though, so this might need to be modified.
                # Also, the WPS files are stored big endian.

                record_size = struct.unpack(">i", raw_data[record_size_idx : 
                                            record_size_idx + SIZE_BYTES])
                record_start = record_size_idx + SIZE_BYTES
                record_end = record_start + record_size[0]
                record_data = raw_data[record_start : record_end]

                parsed_record = parse_record(record_idx, record_data, iproj,
                                             nx, ny)

                try:
                    field_name = parsed_record["field"].strip()
                except KeyError:
                    pass
                else:
                    if field is not None:
                        if field_name != field:
                            keep_record = False

                try:
                    iproj = parsed_record["iproj"]
                except KeyError:
                    pass

                try:
                    nx = parsed_record["nx"]
                except KeyError:
                    pass

                try:
                    ny = parsed_record["ny"]
                except KeyError:
                    pass

                wps_params.update(parsed_record)

                record_size_idx = record_end + SIZE_BYTES

            if keep_record:
                record_list.append(WPSData(**wps_params))   

            # Repeat for all record slabs
            if record_end + SIZE_BYTES > end_of_file_idx:
                break

    return record_list

Hope this helps.

[lesserwhirls]

On the Unidata side of things, we've decided to go with directly returning an intermediate file from the TDS rather than returning netCDF and converting that to the intermediate file format using python. @lmadaus had mentioned that they were working on a WRF Intermediate File Format writer in Python over at Jupiter, but I'm not sure of the details.