Skip to content

Commit

Permalink
Fix grid encoding (#12)
Browse files Browse the repository at this point in the history
Need to account for differences in the way times are encoded into
GEMPAK files based on whether they are an analysis or forecast.

Fix issue with analysis grids not having the right time written into the
GEMPAK grid. Bumped version.

Fixed issue when user inputs only a grid type to the date_time string
but not a forecast time. Time will now default to zero.

Removed TVM/UTM as they appear to not be supported since an old GEMPAK
version. Fixed issues with LEA/GNO projections as pyproj (or PROJ?) does
not seem to convert to CF format. Added several additional tests for
reading and writing all grid projections.
  • Loading branch information
nawendt committed Jul 29, 2024
1 parent c415f30 commit 15486b7
Show file tree
Hide file tree
Showing 23 changed files with 621 additions and 108 deletions.
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -106,3 +106,6 @@ venv.bak/
### VisualStudioCode ###
.vscode/*
**/.history

### GEMPAK ###
*.nts
46 changes: 21 additions & 25 deletions examples/gempakio_examples.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -558,15 +558,23 @@
"id": "d0bbe9fc",
"metadata": {},
"source": [
"Each grid file will contain data on a single, common grid. Data will often be several atmospheric quantities (moisture, temperature, etc.) that can be at different levels, in different vertical coordinates, and at diferent times. We can make a grid using the dewpoint data that was loaded earlier."
]
},
{
"cell_type": "markdown",
"id": "0303490b",
"metadata": {},
"source": [
"The initial steps to creating a grid will require the coordinates (lat/lon or projected x/y) and the projection information. Unprojected lat/lon should use the equidistant cylindrical projection (CED in GEMPAK). Other possible projections are described in the [GEMPAK documentation](https://unidata.github.io/gempak/man/parm/proj.html). The dimension order for the input coordinates and data will be (y, x)."
"Each grid file will contain data on a single, common grid. Data will often be several atmospheric quantities (moisture, temperature, etc.) that can be at different levels, in different vertical coordinates, and at diferent times. We can make a grid using the dewpoint data that was loaded earlier.\n",
"\n",
"The initial steps to creating a grid will require the coordinates (lat/lon or projected x/y) and the projection information. Unprojected lat/lon should use the equidistant cylindrical projection (CED in GEMPAK). Other possible projections are described in the [GEMPAK documentation](https://unidata.github.io/gempak/man/parm/proj.html). The dimension order for the input coordinates and data will be (y, x). For coordinate and data input, the origin is assumed to be the lower left corner. You may have to transform your grids to fit this convention.\n",
"\n",
"When inputing time(s) for grids, there are different types that can be used. Setting the grid type explicity is only an option when using strings for input (with the format of YYYYmmddHHMMthhhmm). The `t` in the previous date/time string can be set to any of the following characters:\n",
"\n",
"| Grid Type | Code |\n",
"|:--------------|:--------------:|\n",
"|Analysis | A |\n",
"|Forecast | F |\n",
"|Valid | V |\n",
"|Guess | G |\n",
"|Initial | I |\n",
"\n",
"If a `datetime` object is used, an analsyis grid is assumed. The same goes for a date/time string that does not include any of the `thhhmm` portion. When setting the grid type, you'll want to provide, at a minimum, the `hhh` or hours beyond the initial time. One to three digits will work. If you set the grid type but omit the hours, 0 will be assumed. If you want to also set the minutes, you will have to use at least four digits. Being able to set the grid type becomes important when using grids in NMAP. Issues can arise if the grid type or forecast times do not match what is expected based on the file template you use in the `datatype.tbl` configuration.\n",
"\n",
"As an example, if a grid represented a 3-hour forecast from the initialization time of 2012-04-14 18:00 UTC the appropriate date/time string would be `201204141800F003`."
]
},
{
Expand Down Expand Up @@ -628,14 +636,8 @@
"id": "4197c642",
"metadata": {},
"source": [
"Sounding files will contain a set of parameters that all vertical profiles will have. The vertical coordinates used in these file are generally going to be either height or pressure. While the parameters will not change after initialization, the location of the profile and the levels at which the data are can all vary."
]
},
{
"cell_type": "markdown",
"id": "01fdbe35",
"metadata": {},
"source": [
"Sounding files will contain a set of parameters that all vertical profiles will have. The vertical coordinates used in these file are generally going to be either height or pressure. While the parameters will not change after initialization, the location of the profile and the levels at which the data are can all vary.\n",
"\n",
"The parameters that are going to be needed for display in GEMPAK are PRES, HGHT, TEMP, DWPT, DRCT (wind direction), and SPED (wind speed). Other parameters can be added as desired. For more information on sounding parameters, see the [GEMPAK documentation](https://unidata.github.io/gempak/man/parm/snparm.html)."
]
},
Expand Down Expand Up @@ -716,14 +718,8 @@
"id": "1d2220fa",
"metadata": {},
"source": [
"Surface files are similar to sounding files in how they are structured. A set of parameters will be present for each station that is in the file. Stations can be any location and have different times associate with them."
]
},
{
"cell_type": "markdown",
"id": "eafe15ef",
"metadata": {},
"source": [
"Surface files are similar to sounding files in how they are structured. A set of parameters will be present for each station that is in the file. Stations can be any location and have different times associate with them.\n",
"\n",
"Ideally, to display surface data in GEMPAK, you will want to have parameters that correspond to temperature, dewpoint, wind speed, wind direction, and sky cover. The surface data loaded earlier can be used to create a new file. With some clever looping we can add multiple stations. While many GEMPAK surface files have data packed in order to compress and save space, this is not currently implemented in gempakIO. For more information on surface parameters, see the [GEMPAK documentation](https://unidata.github.io/gempak/man/parm/sfparm.html)."
]
},
Expand Down
2 changes: 1 addition & 1 deletion src/gempakio/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,4 +7,4 @@
from gempakio.decode.vgf import VectorGraphicFile
from gempakio.encode.gempak import GridFile, SoundingFile, SurfaceFile

__version__ = '1.1.0'
__version__ = '1.1.1'
32 changes: 13 additions & 19 deletions src/gempakio/decode/gempak.py
Original file line number Diff line number Diff line change
Expand Up @@ -49,9 +49,7 @@
'AED': ('aeqd', 'azm'),
'ORT': ('ortho', 'azm'),
'LEA': ('laea', 'azm'),
'GNO': ('gnom', 'azm'),
'TVM': ('tmerc', 'obq'),
'UTM': ('utm', 'obq'),
'GNO': ('gnom', 'azm')
}
GVCORD_TO_VAR = {
'PRES': 'p',
Expand Down Expand Up @@ -341,7 +339,12 @@ def _process_gempak_header(self):

@staticmethod
def _convert_dattim(dattim):
"""Convert GEMPAK DATTIM integer to datetime object."""
"""Convert GEMPAK DATTIM integer to datetime object.
Notes
-----
See GEMPAK subroutine TG_FTOI.
"""
if dattim:
if dattim < 100000000:
dt = datetime.strptime(f'{dattim:06d}', '%y%m%d')
Expand All @@ -353,7 +356,12 @@ def _convert_dattim(dattim):

@staticmethod
def _convert_ftime(ftime):
"""Convert GEMPAK forecast time and type integer."""
"""Convert GEMPAK forecast time and type integer.
Notes
-----
See GEMPAK subroutine TG_CFTM.
"""
if ftime >= 0:
iftype = ForecastType(ftime // 100000)
iftime = ftime - iftype.value * 100000
Expand Down Expand Up @@ -609,20 +617,6 @@ def _get_crs(self):
'ellps': ellps,
'R': earth_radius})

elif ptype == 'obq':
lon_0 = self.navigation_block.proj_angle1
if gemproj == 'UTM':
zone = np.digitize((lon_0 % 360) / 6 + 1, range(1, 61), right=True)
self.crs = pyproj.CRS.from_dict({'proj': proj,
'zone': zone,
'ellps': ellps,
'R': earth_radius})
else:
self.crs = pyproj.CRS.from_dict({'proj': proj,
'lon_0': lon_0,
'ellps': ellps,
'R': earth_radius})

def _set_coordinates(self):
"""Use GEMPAK navigation block to define coordinates.
Expand Down
Loading

0 comments on commit 15486b7

Please sign in to comment.