{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Read Fortran Binary Data Files"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Although most of the data files I work with in atmosphere, ocean, and climate science are in self describing formats, such as netCDF, grib, etc, I still sometimes encounter binary data files written in fortran direct or sequential access. I want to be able to read in these files and convert them to `xarray.Dataset` so I can use all the `xarray` methods and tools to process them.\n",
"\n",
"This notebook will demonstrate how to read a fortran binary sequential access file and convert it to an `xarray.Dataset`\n",
"\n",
"Thanks to Phil Pegion at University of Colorado/CIRES and NOAA/PSL for showing me how to do this."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I was given the location of the following file from a colleague (Thanks David Straus!). It consists of EOFs (spatial patterns of variability) based on anomalies of 500 hPa geopotential height and 250 hPa zonal winds. The file is located here on the COLA servers:\n",
"```/project/mjo/straus/ERAI_T42/eofs/DJF/5day_means/eofs_Z500U250_PNA_DJF_T42.1980-2014.dat```\n",
"\n",
"A [GrADS](http://cola.gmu.edu/grads/) .ctl file was provided with the data. It has the important metadata that will be needed to know how to read the file. Here are its contents:\n",
"\n",
"```\n",
"dset /project/mjo/straus/ERAI_T42/eofs/DJF/5day_means/eofs_Z500U250_PNA_DJF_T42.1980-2014.dat\n",
"undef -9999.9\n",
"title T42 gridded Minerva EOFs config Z500U250_PNA_DJF_T42\n",
"options sequential yrev \n",
"xdef 128 linear 0.0 2.8125\n",
"ydef 64 LEVELS -87.86 -85.10 -82.31 -79.53 -76.74 -73.95 -71.16 -68.37\n",
" -65.58 -62.79 -60.00 -57.21 -54.42 -51.63 -48.84 -46.04 -43.25 -40.46\n",
" -37.67 -34.88 -32.09 -29.30 -26.51 -23.72 -20.93 -18.14 -15.35 -12.56\n",
" -9.77 -6.98 -4.19 -1.40 1.40 4.19 6.98 9.77 12.56 15.35\n",
" 18.14 20.93 23.72 26.51 29.30 32.09 34.88 37.67 40.46 43.25\n",
" 46.04 48.84 51.63 54.42 57.21 60.00 62.79 65.58 68.37 71.16\n",
" 73.95 76.74 79.53 82.31 85.10 87.86\n",
"zdef 1 levels 1000\n",
"tdef 50 linear 01dec1980 1dy\n",
"vars 2\n",
"geo 0 9 geo 500\n",
"uwn 0 9 geo 250\n",
"endvars\n",
"```\n",
"\n",
"The key things to note from the .ctl file are:\n",
"* lons = 128, start from 0.0 and increase in increments of 2.8125\n",
"* lats = 64 and are specified\n",
"* the file has 50 \"times\"; these correspond to 50 EOFs\n",
"* there are 2 variables in the file, geo (500 hPa geopotential height) and uwn (250 hPa zonal wind)\n",
"* under options, it says sequential, so our file is sequential access. This means that each record has 2 extra INTEGER*4 in it. One at the beginning of the record and one at the end of the record.\n",
"* Missing data is indicated with values of -9999.9\n",
"\n",
"The format of the data file is based on the format that GrADS expects, described [here (http://cola.gmu.edu/grads/gadoc/aboutgriddeddata.html#structure). \n",
"\n",
"Each record is a grid of all lats and lons. The records are in the following order:\n",
"\n",
"1. Time 1, geo\n",
"2. Time 1, uwn\n",
"3. Time 2, geo\n",
"4. Time 2, uwn\n",
"5. ..."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from array import array\n",
"import xarray as xr\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the path and filename"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"path='/project/mjo/straus/ERAI_T42/eofs/DJF/5day_means/'\n",
"fname='eofs_Z500U250_PNA_DJF_T42.1980-2014.dat'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the dimensions of the data based on the .ctl file"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [],
"source": [
"nlons=128\n",
"nlats=64\n",
"neofs=50\n",
"nvars=2\n",
"missing_value=-9999.9"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the coordinates of the data based on the .ctl file"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [],
"source": [
"# Define lats\n",
"lats=[-87.86, -85.10, -82.31, -79.53, -76.74, -73.95, -71.16, -68.37, \\\n",
" -65.58, -62.79, -60.00, -57.21, -54.42, -51.63, -48.84, -46.04, \\\n",
" -43.25, -40.46, -37.67, -34.88, -32.09, -29.30, -26.51, -23.72, \\\n",
" -20.93, -18.14, -15.35, -12.56, -9.77, -6.98, -4.19, -1.40, \\\n",
" 1.40, 4.19, 6.98, 9.77, 12.56, 15.35, 18.14, 20.93, \\\n",
" 23.72, 26.51, 29.30, 32.09, 34.88, 37.67, 40.46, 43.25, \\\n",
" 46.04, 48.84, 51.63, 54.42, 57.21, 60.00, 62.79, 65.58, \\\n",
" 68.37, 71.16, 73.95, 76.74, 79.53, 82.31, 85.10, 87.86]\n",
"\n",
"# Yrev in .ctl file indicates lats are reversed\n",
"lats=lats[::-1]\n",
"\n",
"# Define lons\n",
"lons=np.arange(128)*2.8125 + 0.0\n",
"\n",
"# Define times as a pandas date range\n",
"eofs=np.arange(neofs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the length of each record. Be sure to include the 2 extra integers for sequential access."
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"recl=(nlons*nlats+2)*4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create empty array to store the data"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"data=np.zeros((neofs,nlats,nlons,nvars))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Read the data"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"# Open file\n",
"luin = open(path+fname,'rb')\n",
"\n",
"# Loop over all times\n",
"for e in range(neofs):\n",
" \n",
" # Loop over both variables\n",
" for v in range(nvars):\n",
" \n",
" # Read in fortran record in bytes\n",
" tmp=luin.read(recl) \n",
" \n",
" # Convert to single precision (real 32bit)\n",
" tmp1=array('f',tmp) \n",
" \n",
" # Pull out data array (leaving behind fortran control records)for fortran sequential\n",
" tmp2=tmp1[1:-1] \n",
" \n",
" # Create a 2d array (lat x lon) and store it in the data array\n",
" data[e,:,:,v]=np.reshape(tmp2,(nlats,nlons)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Extract out our two variables"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [],
"source": [
"z500=data[:,:,:,0]\n",
"u250=data[:,:,:,1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Take care of missing data by setting it to NAN"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"z500[z500<=missing_value]=np.nan\n",
"u250[u250<=missing_value]=np.nan"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Put the data into an `xarray.Dataset`"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"# 500 hPa Geopotential Height\n",
"z500_ds=xr.DataArray(z500,\n",
" coords={'eofnum':eofs,\n",
" 'lat':lats,\n",
" 'lon': lons},\n",
" dims=['eofnum','lat','lon']) \n",
"z500_ds=z500_ds.to_dataset(name='z500')\n",
"\n",
"# 250 hPa Zonal Wind\n",
"u250_ds=xr.DataArray(u250,\n",
" coords={'eofnum':eofs,\n",
" 'lat':lats,\n",
" 'lon': lons},\n",
" dims=['eofnum','lat','lon']) \n",
"u250_ds=u250_ds.to_dataset(name='u250')\n",
"\n",
"# Merge to have both in the same `xarray.Dataset`\n",
"ds=xr.merge([z500_ds,u250_ds])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This dataset has global values, but only contains valid data in a certain region. All other values are marked as missing. So, we will drop missing data for `lat` and `lon` where `all` the data are missing."
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [],
"source": [
"ds=ds.dropna(dim='lon',how='all').dropna(dim='lat',how='all')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have an `xarray.Dataset` with both variables in it"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"