<jupyter_start><jupyter_text>netCDF File Visualization Case StudyI was asked by a colleague to visualize data contained within this [netCDF file](https://motherlode.ucar.edu/repository/entry/show/RAMADDA/Unidata/Staff/Julien+Chastang/netcdf-explore?entryid=c7239224-d3fe-45d8-b100-43ae043824c3) ([OPeNDAP link](https://motherlode.ucar.edu/repository/opendap/41f2b38a-4e70-4135-8ff8-dbf3d1dcbfc1/entry.das)) with Python. What follows is an exploration of how I achieved that objective. Because this exercise touches upon many technologies related to Unidata, it makes for an interesting case study. We will be meandering through,- netCDF- WMO GRIB metadata- Map projections- xarray data analysis library- cartopy visualization library  Crack Open the FileTo get our bearings let's see what there is inside our netCDF file. We will be using the [xarray library](https://github.com/pydata/xarray) to dig inside our netCDF data. xarray is similar to pandas, but for the [Common Data Model](http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/CDM/). We could have just used the [netcdf4-python library](https://github.com/Unidata/netcdf4-python) but xarray has output that is more nicely formatted. Let's first import xarray and open the dataset.<jupyter_code>import xarray
ds = xarray.open_dataset('https://motherlode.ucar.edu/repository/opendap/41f2b38a-4e70-4135-8ff8-dbf3d1dcbfc1/entry.das', 
                       decode_times=False)
print(ds)<jupyter_output><xarray.Dataset>
Dimensions:         (nav: 1, record: 1, x: 1901, y: 1801, z: 1)
Coordinates:
  * z               (z) float32 0.0
Dimensions without coordinates: nav, record, x, y
Data variables:
    u               (record, z, y, x) float32 ...
    v               (record, z, y, x) float32 ...
    p               (record, z, y, x) float32 ...
    t               (record, z, y, x) float32 ...
    td              (record, z, y, x) float32 ...
    vv              (record, z, y, x) float32 ...
    rh              (record, z, y, x) float32 ...
    msl             (record, z, y, x) float32 ...
    tad             (record, z, y, x) float32 ...
    th              (record, z, y, x) float32 ...
    the             (record, z, y, x) float32 ...
    ps              (record, z, y, x) float32 ...
    vor             (record, z, y, x) float32 ...
    mr              (record, z, y, x) float32 ...
    mrc             (record, z, y, x) float32 ...
    div             (record, z, y, x) float32 ...
    [...]<jupyter_text>Dimensions, Coordinates, Data VariablesAs far as the dimensions and coordinates go, the most relevant and important coordinates variables are `x` and `y`. We can see the data variables, such as, temperature (`t`), mixing ratio (`mr`), and potential temperature (`th`), are mostly on a 1901 x 1801 grid. There is also the mysterious `nav` dimension and associated data variables which we will be examining later.Let's set a goal of visualizing **potential temperature** with the [Cartopy](http://scitools.org.uk/cartopy/) plotting package.The first step is to get more information concerning the variables we are interested in. For example, let's look at _potential temperature_ or `th`.<jupyter_code>print(ds['th'])<jupyter_output><xarray.DataArray 'th' (record: 1, z: 1, y: 1801, x: 1901)>
[3423701 values with dtype=float32]
Coordinates:
  * z        (z) float32 0.0
Dimensions without coordinates: record, y, x
Attributes:
    navigation_dim:  nav
    record:          valtime, reftime
    long_name:       potential temperature
    units:           kelvin
    valid_range:     [-75. 125.]
    LAPS_var:        TH
    lvl_coord:       AGL
    LAPS_units:      K<jupyter_text>potential temperature (`th`)Let's grab the data array for potential temperature (`th`).<jupyter_code>th = ds['th'].values[0][0]
print(th)<jupyter_output>[[282.0028  281.97714 281.89117 ... 283.29944 281.8932  281.56363]
 [282.03864 282.0177  281.9361  ... 282.73636 281.5668  281.56363]
 [282.05988 282.05298 282.01532 ... 282.0037  281.56363 281.5663 ]
 ...
 [283.70425 283.44934 283.6949  ... 281.56363 281.56363 281.56363]
 [283.87967 283.49844 283.4708  ... 281.56363 281.56363 281.56363]
 [283.983   284.04346 284.17642 ... 281.56363 281.56363 281.56363]]<jupyter_text>To Visualize the Data, We have to Decrypt the ProjectionIn order, to visualize the data that are contained within a two-dimensional array onto a map that represents a three-dimensional globe, we need to understand the projection of the data.We can make an educated guess these are contained in the data variables with the `nav` cooridinate variable.Specifically,- `grid_type`- `grid_type_code`- `x_dim`- `y_dim`- `Nx`- `Ny`- `La1`- `Lo1`- `LoV`- `Latin1`- `Latin2`- `Dx`- `Dy`**But what are these??**  For Grins, Let's Scrutinize the `grid_type_code`<jupyter_code>print(ds['grid_type_code'])<jupyter_output><xarray.DataArray 'grid_type_code' (nav: 1)>
array([5], dtype=int32)
Dimensions without coordinates: nav
Attributes:
    long_name:  GRIB-1 GDS data representation type<jupyter_text>Google to the RescueA simple Google search of `GRIB-1 GDS data representation type` takes us to[A GUIDE TO THE CODE FORM FM 92-IX Ext. GRIB Edition 1 from 1994](http://www.wmo.int/pages/prog/www/WMOCodes/Guides/GRIB/GRIB1-Contents.html "GRIB") document. Therein one can find an explanation of the variables needed to understand the map projection. Let's review these variables.<jupyter_code>print(ds['grid_type_code'].values[0])<jupyter_output>5<jupyter_text>What is `grid_type_code` of `5`?Let's look at [Table 6 ](http://www.wmo.int/pages/prog/www/WMOCodes/Guides/GRIB/GRIB1-Contents.html "GRIB Projection Definitions"). A `grid_type_code` of `5` corresponds to a projection of **Polar Stereographic**.  Next up `grid_type`<jupyter_code>grid_type = ds['grid_type'].values
print('The grid type is ', grid_type[0])<jupyter_output>The grid type is  b'secant lambert conformal      '<jupyter_text>Uh oh! Polar Stereographic or Lambert Conformal??_Note that this newest piece of information relating to a Lambert Conformal projection disagrees with the  earlier projection information about a Polar Stereographic projection._ There is a **bug** in the metadata description of the projection.  Moving on Anyway, next `Nx` and `Ny`According to the grib documentation `Nx` and `Ny` represent the number grid points along the x and y axes. Let's grab those.<jupyter_code>nx, ny = ds['Nx'].values[0], ds['Ny'].values[0]
print(nx, ny)<jupyter_output>1901 1801<jupyter_text>`La1` and `Lo1` Next let's get `La1` and `Lo1` which are defined as the "first grid points" These are probably the latitude and longitude for one of the corners of the grid.<jupyter_code>la1, lo1 = ds['La1'].values[0], ds['Lo1'].values[0]
print(la1, lo1)<jupyter_output>10.160843 78.986084<jupyter_text>`Latin1` and `Latin2`Next up are the rather mysteriously named `Latin1` and `Latin2` variables. When I first saw these identifiers, I thought they referred to a Unicode block, but in fact they relate to the secants of the projection cone. I do not know why they are called "Latin" and this name is confusing. **At any rate, we can feel comfortable that we are dealing with Lambert Conformal rather than Polar Stereographic.**![Lambert Conformal](http://www.geo.hunter.cuny.edu/~jochen/gtech201/Lectures/Lec6concepts/Map%20coordinate%20systems/Lambert%20Conformal%20Conic_files/image002.gif "Lambert Conformal")Credit: http://www.geo.hunter.cuny.edu/~jochen<jupyter_code>latin1, latin2 = ds['Latin1'].values[0], ds['Latin2'].values[0]
print(latin1, latin2)<jupyter_output>30.0 60.0<jupyter_text>The Central Meridian for the Lambert Conformal Projection, `LoV`If we are defining a Lambert Conformal projection, we will require the central meridian that the GRIB documentation refers to as `LoV`.<jupyter_code>lov = ds['LoV'].values[0]
print(lov)<jupyter_output>102.0<jupyter_text>`Dx` and `Dy`Finally, let's look at the grid increments. In particular, we need to find the units.<jupyter_code>print(ds['Dx'])
print(ds['Dy'])<jupyter_output><xarray.DataArray 'Dx' (nav: 1)>
array([3000.], dtype=float32)
Dimensions without coordinates: nav
Attributes:
    long_name:  x grid increment
    units:      meters
<xarray.DataArray 'Dy' (nav: 1)>
array([3000.], dtype=float32)
Dimensions without coordinates: nav
Attributes:
    long_name:  y grid increment
    units:      meters<jupyter_text>Units for `Dx` and `Dy`The units for the deltas are in meters.<jupyter_code>dx,dy = ds['Dx'].values[0],ds['Dy'].values[0]
print(dx,dy)<jupyter_output>3000.0 3000.0<jupyter_text>Let's Review What We HaveWe now have all the information we need to understand the Lambert projection:- The secants of the Lambert Conformal projection (`Latin1`, `Latin2`)- The central meridian of the projection (`LoV`)Moreover, we have additional information that shows how the data grid relates to the projection:    - The number of grid points in x and y (`Nx`, `Ny`)- The delta in meters between grid point (`Dx`, `Dy`)- The first latitude and longitude of the data (`first latitude`, `first longitude`).  We are Ready for Visualization (almost)!Let's import **cartopy** and **matplotlib**.<jupyter_code>%matplotlib inline

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import matplotlib as mpl<jupyter_output><empty_output><jupyter_text>Define the Lambert Conformal Projection with Cartopy<jupyter_code>proj = ccrs.LambertConformal(central_longitude=lov,standard_parallels=(latin1,latin2))<jupyter_output><empty_output><jupyter_text>Lambert Conformal Grid Extents- To plot the data we need the `left`,`right`,`bottom`,`top` extents of the grid **expressed in Lambert Conformalcoordinates**.- __**Key point**: The projection coordinate systems have flat topology and Euclidean distance.__   Calculating the ExtentsRemember, we have:- The number of grid points in x and y (`Nx`, `Ny`)- The delta in meters between grid point (`Dx`, `Dy`)- The first latitude and longitude of the data (`first latitude`, `first longitude`).We have one of the corners in latitude and longitude, but we need to convert it LC coordinates and derive the other corner.  Platte Carrée ProjectionThe Platte Carrée Projection is a very simple X,Y/Cartesian projection. It is used a lot in Cartopy because it allows you to express coordinates in familiar Latitude and Longitudes. **Remember**: The projection coordinate systems have flat topology and Euclidean distance.  Platte Carrée![Platte Carree](https://upload.wikimedia.org/wikipedia/commons/8/83/Equirectangular_projection_SW.jpg "Platte Carree")Source: [Wikipedia Source](https://en.wikipedia.org/wiki/Equirectangular_projection)  Create the PlatteCarre Cartopy Projection<jupyter_code>pc = ccrs.PlateCarree()<jupyter_output><empty_output><jupyter_text>Convert Corner from Lat/Lon PlatteCarre to LCThe `transform_point` method translates coordinates from one projection coordinate system to the other.<jupyter_code>left,bottom = proj.transform_point(lo1,la1,pc)
print(left,bottom)<jupyter_output>-2851705.4600501494 -2903564.859896133<jupyter_text>Derive Opposite CornerDerive the opposite corner from the number of points and the delta. **Again**, we can do this because the projection coordinate systems have flat topology and Euclidean distance.<jupyter_code>right,top = left + nx*dx,bottom + ny*dy
print(right,top)<jupyter_output>2851294.5399498506 2499435.140103867<jupyter_text>Plot It Up!We now have the extents, we are ready to plot.<jupyter_code>#Define the figure
fig = plt.figure(figsize=(12, 12))

# Define the extents and add the data
ax = plt.axes(projection=proj)
extents = (left, right, bottom, top)
ax.contourf(th, origin='lower', extent=extents, transform=proj)

# Add bells and whistles
ax.coastlines(resolution='50m', color='black', linewidth=2)
ax.add_feature(ccrs.cartopy.feature.STATES)
ax.add_feature(ccrs.cartopy.feature.BORDERS, linewidth=1, edgecolor='black')
ax.gridlines()

plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1087593072868736, 'fraction_numerical': 0.05279009453693633, 'mean_word_length': 4.465935070873343, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 64, 'https://': 7, 'lorem ipsum': 0, 'www.': 7, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '256514', 'license': 'MIT', 'n_tokens_mistral': 4344, 'n_tokens_neox': 3803, 'n_words': 1271, 'path': 'failing_notebooks/casestudy.ipynb', 'repo_name': 'fuyihang98/python-workshop', 'size': '11953'}