{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "(nb_testing_metpy)=\n", "# Getting started with MetPy\n", "\n", "Before we get started, we test the learning environment and the most important packages needed to run the notebooks. This is not so much a continuous coherent exercise as individual examples based on the different packages.This exercise is neither an introduction to Python nor extensive tutorials for the individual packages. I advise you, if you have little or no experience with the packages, to work through the relevant tutorial on the websites. All packages offer very good and extensive tutorials. Most of the functions presented here have been taken from these websites." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Learning objectives:
\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "
\n", "How to proceed:
\n", " \n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MetPy is a collection of tools in Python for reading, visualizing, and performing calculations with weather data. One of the most significant differences in syntax for MetPy, compared to other Python libraries, is the frequent requirement of units to be attached to arrays before being passed to MetPy functions. There are very few exceptions to this, and you’ll usually be safer to always use units whenever applicable to make sure that your analyses are done correctly. Once you get used to the units syntax, it becomes very handy, as you never have to worry about unit conversion for any calculation. MetPy does it for you!\n", "\n", "Let's load the MetPy and numpy library\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Load the pandas package\n", "import numpy as np\n", "from metpy.units import units" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Please note, we have loaded only the units module here. We can assign a unit to an array by" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Attach the unit meters to the distance array\n", "distance = np.arange(1, 5) * units.meters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we can attach create a time array" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# This is another way to attach units to an array \n", "time = units.Quantity(np.arange(2, 10, 2), 'sec')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can simply do unit-aware calculations" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Magnitude
[0.5 0.5 0.5 0.5]
Unitsmeter/second
" ], "text/latex": [ "$\\begin{pmatrix}0.5 & 0.5 & 0.5 & 0.5\\end{pmatrix}\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$" ], "text/plain": [ "array([0.5, 0.5, 0.5, 0.5]) " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Calculate the velocity\n", "distance/time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## MetPy with xarray\n", "\n", "MetPy works great with xarray. MetPy’s suite of meteorological calculations are designed to integrate with xarray DataArrays and provides DataArray and Dataset accessors (collections of methods and properties attached to the .metpy property) for coordinate/CRS and unit operations.\n", "\n", "First, some imports ..." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Import xarray\n", "import xarray as xr\n", "\n", "# Any import of metpy activates the accessors\n", "import metpy\n", "from metpy.cbook import get_test_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can open a netCDF-file using xarray" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:                              (time1: 9, latitude: 81,\n",
       "                                          isobaric3: 31, isobaric1: 21,\n",
       "                                          longitude: 131)\n",
       "Coordinates:\n",
       "  * time1                                (time1) datetime64[ns] 2017-09-05T12...\n",
       "    reftime                              datetime64[ns] ...\n",
       "  * latitude                             (latitude) float32 50.0 49.5 ... 10.0\n",
       "  * isobaric3                            (isobaric3) float64 100.0 ... 1e+05\n",
       "  * isobaric1                            (isobaric1) float64 1e+04 ... 1e+05\n",
       "  * longitude                            (longitude) float32 250.0 ... 315.0\n",
       "Data variables:\n",
       "    Vertical_velocity_pressure_isobaric  (time1, isobaric1, latitude, longitude) float32 ...\n",
       "    Relative_humidity_isobaric           (time1, isobaric3, latitude, longitude) float32 ...\n",
       "    Temperature_isobaric                 (time1, isobaric3, latitude, longitude) float32 ...\n",
       "    u-component_of_wind_isobaric         (time1, isobaric3, latitude, longitude) float32 ...\n",
       "    v-component_of_wind_isobaric         (time1, isobaric3, latitude, longitude) float32 ...\n",
       "    Geopotential_height_isobaric         (time1, isobaric3, latitude, longitude) float32 ...\n",
       "    LatLon_361X720-0p25S-180p00E         int32 ...\n",
       "Attributes: (12/13)\n",
       "    Originating_or_generating_Center:                                        ...\n",
       "    Originating_or_generating_Subcenter:                                     ...\n",
       "    GRIB_table_version:                                                      ...\n",
       "    Type_of_generating_process:                                              ...\n",
       "    Analysis_or_forecast_generating_process_identifier_defined_by_originating...\n",
       "    Conventions:                                                             ...\n",
       "    ...                                                                                ...\n",
       "    featureType:                                                             ...\n",
       "    History:                                                                 ...\n",
       "    geospatial_lat_min:                                                      ...\n",
       "    geospatial_lat_max:                                                      ...\n",
       "    geospatial_lon_min:                                                      ...\n",
       "    geospatial_lon_max:                                                      ...
" ], "text/plain": [ "\n", "Dimensions: (time1: 9, latitude: 81,\n", " isobaric3: 31, isobaric1: 21,\n", " longitude: 131)\n", "Coordinates:\n", " * time1 (time1) datetime64[ns] 2017-09-05T12...\n", " reftime datetime64[ns] ...\n", " * latitude (latitude) float32 50.0 49.5 ... 10.0\n", " * isobaric3 (isobaric3) float64 100.0 ... 1e+05\n", " * isobaric1 (isobaric1) float64 1e+04 ... 1e+05\n", " * longitude (longitude) float32 250.0 ... 315.0\n", "Data variables:\n", " Vertical_velocity_pressure_isobaric (time1, isobaric1, latitude, longitude) float32 ...\n", " Relative_humidity_isobaric (time1, isobaric3, latitude, longitude) float32 ...\n", " Temperature_isobaric (time1, isobaric3, latitude, longitude) float32 ...\n", " u-component_of_wind_isobaric (time1, isobaric3, latitude, longitude) float32 ...\n", " v-component_of_wind_isobaric (time1, isobaric3, latitude, longitude) float32 ...\n", " Geopotential_height_isobaric (time1, isobaric3, latitude, longitude) float32 ...\n", " LatLon_361X720-0p25S-180p00E int32 ...\n", "Attributes: (12/13)\n", " Originating_or_generating_Center: ...\n", " Originating_or_generating_Subcenter: ...\n", " GRIB_table_version: ...\n", " Type_of_generating_process: ...\n", " Analysis_or_forecast_generating_process_identifier_defined_by_originating...\n", " Conventions: ...\n", " ... ...\n", " featureType: ...\n", " History: ...\n", " geospatial_lat_min: ...\n", " geospatial_lat_max: ...\n", " geospatial_lon_min: ...\n", " geospatial_lon_max: ..." ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Open the netCDF file as a xarray Dataset\n", "ds = xr.open_dataset(get_test_data('irma_gfs_example.nc', as_file_obj = False))\n", "\n", "# View a summary of the Dataset\n", "ds" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This Dataset consists of dimensions and their associated coordinates, which in turn make up the axes along which the data variables are defined. The dataset also has a dictionary-like collection of attributes. What happens if we look at just a single data variable?" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'Temperature_isobaric' (time1: 9, isobaric3: 31,\n",
       "                                          latitude: 81, longitude: 131)>\n",
       "[2960469 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * time1      (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...\n",
       "    reftime    datetime64[ns] ...\n",
       "  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n",
       "  * isobaric3  (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05\n",
       "  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n",
       "Attributes:\n",
       "    long_name:                      Temperature @ Isobaric surface\n",
       "    units:                          K\n",
       "    Grib_Variable_Id:               VAR_0-0-0_L100\n",
       "    Grib2_Parameter:                [0 0 0]\n",
       "    Grib2_Parameter_Discipline:     Meteorological products\n",
       "    Grib2_Parameter_Category:       Temperature\n",
       "    Grib2_Parameter_Name:           Temperature\n",
       "    Grib2_Level_Type:               100\n",
       "    Grib2_Level_Desc:               Isobaric surface\n",
       "    Grib2_Generating_Process_Type:  Forecast\n",
       "    grid_mapping:                   LatLon_361X720-0p25S-180p00E
" ], "text/plain": [ "\n", "[2960469 values with dtype=float32]\n", "Coordinates:\n", " * time1 (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...\n", " reftime datetime64[ns] ...\n", " * latitude (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n", " * isobaric3 (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05\n", " * longitude (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n", "Attributes:\n", " long_name: Temperature @ Isobaric surface\n", " units: K\n", " Grib_Variable_Id: VAR_0-0-0_L100\n", " Grib2_Parameter: [0 0 0]\n", " Grib2_Parameter_Discipline: Meteorological products\n", " Grib2_Parameter_Category: Temperature\n", " Grib2_Parameter_Name: Temperature\n", " Grib2_Level_Type: 100\n", " Grib2_Level_Desc: Isobaric surface\n", " Grib2_Generating_Process_Type: Forecast\n", " grid_mapping: LatLon_361X720-0p25S-180p00E" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temperature = ds['Temperature_isobaric']\n", "temperature" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "This is a DataArray, which stores just a single data variable with its associated coordinates and attributes. These individual DataArrays are the kinds of objects that MetPy’s calculations take as input (more on that in Calculations section below)." ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Coordinates and Corrdinate Reference Systems" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "MetPy’s first set of helpers comes with identifying coordinate types. In a given dataset, coordinates can have a variety of different names (time, isobaric1 ...). Following CF conventions, as well as using some fall-back regular expressions, MetPy can systematically identify coordinates of the following types:\n", "\n", "
    \n", "
  • time\n", "
  • vertical\n", "
  • latitude\n", "
  • y\n", "
  • longitude\n", "
  • x\n", "
\n", "\n", "When identifying a single coordinate, it is best to use the property directly associated with that type\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'time1' (time1: 9)>\n",
       "array(['2017-09-05T12:00:00.000000000', '2017-09-05T15:00:00.000000000',\n",
       "       '2017-09-05T18:00:00.000000000', '2017-09-05T21:00:00.000000000',\n",
       "       '2017-09-06T00:00:00.000000000', '2017-09-06T03:00:00.000000000',\n",
       "       '2017-09-06T06:00:00.000000000', '2017-09-06T09:00:00.000000000',\n",
       "       '2017-09-06T12:00:00.000000000'], dtype='datetime64[ns]')\n",
       "Coordinates:\n",
       "  * time1    (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:00:00\n",
       "    reftime  datetime64[ns] ...\n",
       "Attributes:\n",
       "    standard_name:  time\n",
       "    long_name:      time\n",
       "    udunits:        Hour since 2017-09-05T12:00:00Z\n",
       "    _metpy_axis:    time
" ], "text/plain": [ "\n", "array(['2017-09-05T12:00:00.000000000', '2017-09-05T15:00:00.000000000',\n", " '2017-09-05T18:00:00.000000000', '2017-09-05T21:00:00.000000000',\n", " '2017-09-06T00:00:00.000000000', '2017-09-06T03:00:00.000000000',\n", " '2017-09-06T06:00:00.000000000', '2017-09-06T09:00:00.000000000',\n", " '2017-09-06T12:00:00.000000000'], dtype='datetime64[ns]')\n", "Coordinates:\n", " * time1 (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:00:00\n", " reftime datetime64[ns] ...\n", "Attributes:\n", " standard_name: time\n", " long_name: time\n", " udunits: Hour since 2017-09-05T12:00:00Z\n", " _metpy_axis: time" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temperature.metpy.time" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'longitude' (longitude: 131)>\n",
       "array([250. , 250.5, 251. , 251.5, 252. , 252.5, 253. , 253.5, 254. , 254.5,\n",
       "       255. , 255.5, 256. , 256.5, 257. , 257.5, 258. , 258.5, 259. , 259.5,\n",
       "       260. , 260.5, 261. , 261.5, 262. , 262.5, 263. , 263.5, 264. , 264.5,\n",
       "       265. , 265.5, 266. , 266.5, 267. , 267.5, 268. , 268.5, 269. , 269.5,\n",
       "       270. , 270.5, 271. , 271.5, 272. , 272.5, 273. , 273.5, 274. , 274.5,\n",
       "       275. , 275.5, 276. , 276.5, 277. , 277.5, 278. , 278.5, 279. , 279.5,\n",
       "       280. , 280.5, 281. , 281.5, 282. , 282.5, 283. , 283.5, 284. , 284.5,\n",
       "       285. , 285.5, 286. , 286.5, 287. , 287.5, 288. , 288.5, 289. , 289.5,\n",
       "       290. , 290.5, 291. , 291.5, 292. , 292.5, 293. , 293.5, 294. , 294.5,\n",
       "       295. , 295.5, 296. , 296.5, 297. , 297.5, 298. , 298.5, 299. , 299.5,\n",
       "       300. , 300.5, 301. , 301.5, 302. , 302.5, 303. , 303.5, 304. , 304.5,\n",
       "       305. , 305.5, 306. , 306.5, 307. , 307.5, 308. , 308.5, 309. , 309.5,\n",
       "       310. , 310.5, 311. , 311.5, 312. , 312.5, 313. , 313.5, 314. , 314.5,\n",
       "       315. ], dtype=float32)\n",
       "Coordinates:\n",
       "    reftime    datetime64[ns] ...\n",
       "  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n",
       "Attributes:\n",
       "    units:          degrees_east\n",
       "    standard_name:  longitude\n",
       "    _metpy_axis:    x,longitude
" ], "text/plain": [ "\n", "array([250. , 250.5, 251. , 251.5, 252. , 252.5, 253. , 253.5, 254. , 254.5,\n", " 255. , 255.5, 256. , 256.5, 257. , 257.5, 258. , 258.5, 259. , 259.5,\n", " 260. , 260.5, 261. , 261.5, 262. , 262.5, 263. , 263.5, 264. , 264.5,\n", " 265. , 265.5, 266. , 266.5, 267. , 267.5, 268. , 268.5, 269. , 269.5,\n", " 270. , 270.5, 271. , 271.5, 272. , 272.5, 273. , 273.5, 274. , 274.5,\n", " 275. , 275.5, 276. , 276.5, 277. , 277.5, 278. , 278.5, 279. , 279.5,\n", " 280. , 280.5, 281. , 281.5, 282. , 282.5, 283. , 283.5, 284. , 284.5,\n", " 285. , 285.5, 286. , 286.5, 287. , 287.5, 288. , 288.5, 289. , 289.5,\n", " 290. , 290.5, 291. , 291.5, 292. , 292.5, 293. , 293.5, 294. , 294.5,\n", " 295. , 295.5, 296. , 296.5, 297. , 297.5, 298. , 298.5, 299. , 299.5,\n", " 300. , 300.5, 301. , 301.5, 302. , 302.5, 303. , 303.5, 304. , 304.5,\n", " 305. , 305.5, 306. , 306.5, 307. , 307.5, 308. , 308.5, 309. , 309.5,\n", " 310. , 310.5, 311. , 311.5, 312. , 312.5, 313. , 313.5, 314. , 314.5,\n", " 315. ], dtype=float32)\n", "Coordinates:\n", " reftime datetime64[ns] ...\n", " * longitude (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n", "Attributes:\n", " units: degrees_east\n", " standard_name: longitude\n", " _metpy_axis: x,longitude" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, y = temperature.metpy.coordinates('x', 'y')\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These coordinate type aliases can also be used in MetPy’s wrapped .sel and .loc for indexing and selecting on DataArrays. For example, to access 500 hPa heights at 1800Z," ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'Geopotential_height_isobaric' (latitude: 81, longitude: 131)>\n",
       "[10611 values with dtype=float32]\n",
       "Coordinates:\n",
       "    time1      datetime64[ns] 2017-09-05T18:00:00\n",
       "    reftime    datetime64[ns] ...\n",
       "  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n",
       "    isobaric3  float64 5e+04\n",
       "  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n",
       "Attributes:\n",
       "    long_name:                      Geopotential height @ Isobaric surface\n",
       "    units:                          gpm\n",
       "    Grib_Variable_Id:               VAR_0-3-5_L100\n",
       "    Grib2_Parameter:                [0 3 5]\n",
       "    Grib2_Parameter_Discipline:     Meteorological products\n",
       "    Grib2_Parameter_Category:       Mass\n",
       "    Grib2_Parameter_Name:           Geopotential height\n",
       "    Grib2_Level_Type:               100\n",
       "    Grib2_Level_Desc:               Isobaric surface\n",
       "    Grib2_Generating_Process_Type:  Forecast\n",
       "    grid_mapping:                   LatLon_361X720-0p25S-180p00E
" ], "text/plain": [ "\n", "[10611 values with dtype=float32]\n", "Coordinates:\n", " time1 datetime64[ns] 2017-09-05T18:00:00\n", " reftime datetime64[ns] ...\n", " * latitude (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n", " isobaric3 float64 5e+04\n", " * longitude (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n", "Attributes:\n", " long_name: Geopotential height @ Isobaric surface\n", " units: gpm\n", " Grib_Variable_Id: VAR_0-3-5_L100\n", " Grib2_Parameter: [0 3 5]\n", " Grib2_Parameter_Discipline: Meteorological products\n", " Grib2_Parameter_Category: Mass\n", " Grib2_Parameter_Name: Geopotential height\n", " Grib2_Level_Type: 100\n", " Grib2_Level_Desc: Isobaric surface\n", " Grib2_Generating_Process_Type: Forecast\n", " grid_mapping: LatLon_361X720-0p25S-180p00E" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "heights = ds['Geopotential_height_isobaric'].metpy.sel(\n", " time='2017-09-05 18:00',\n", " vertical=50000.\n", ")\n", "heights" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Beyond just the coordinates themselves, a common need for both calculations with and plots of geospatial data is knowing the coordinate reference system (CRS) on which the horizontal spatial coordinates are defined. MetPy follows the CF Conventions for its CRS definitions, which it then caches on the metpy_crs coordinate in order for it to persist through calculations and other array operations. There are two ways to do so in MetPy:\n", "\n", "First, if your dataset is already conforming to the CF Conventions, it will have a grid mapping variable that is associated with the other data variables by the grid_mapping attribute. This is automatically parsed via the .parse_cf() method:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Parse full dataset\n", "data_parsed = ds.metpy.parse_cf()\n", "\n", "# Parse subset of dataset\n", "data_subset = ds.metpy.parse_cf([\n", " 'u-component_of_wind_isobaric',\n", " 'v-component_of_wind_isobaric',\n", " 'Vertical_velocity_pressure_isobaric'\n", "])\n", "\n", "# Parse single variable\n", "relative_humidity = ds.metpy.parse_cf('Relative_humidity_isobaric')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the newly added metpy_crs non-dimension coordinate. Now how can we use this in practice? For individual DataArrays, we can access the cartopy and pyproj objects corresponding to this CRS:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2023-04-19T08:34:15.514125\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.7.0, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", "
<cartopy.crs.PlateCarree object at 0x16dcd7f10>
" ], "text/plain": [ "\n", "Name: unknown\n", "Axis Info [cartesian]:\n", "- E[east]: Easting (unknown)\n", "- N[north]: Northing (unknown)\n", "- h[up]: Ellipsoidal height (metre)\n", "Area of Use:\n", "- undefined\n", "Coordinate Operation:\n", "- name: unknown\n", "- method: Equidistant Cylindrical\n", "Datum: unknown\n", "- Ellipsoid: WGS 84\n", "- Prime Meridian: Greenwich" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Cartopy CRS, useful for plotting\n", "relative_humidity.metpy.cartopy_crs" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "\n", "Name: undefined\n", "Axis Info [ellipsoidal]:\n", "- lon[east]: Longitude (degree)\n", "- lat[north]: Latitude (degree)\n", "Area of Use:\n", "- undefined\n", "Datum: undefined\n", "- Ellipsoid: undefined\n", "- Prime Meridian: Greenwich" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# pyproj CRS, useful for projection transformations and forward/backward azimuth and great\n", "# circle calculations\n", "relative_humidity.metpy.pyproj_crs" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Units\n", "\n", "Since unit-aware calculations are a major part of the MetPy library, unit support is a major part of MetPy’s xarray integration!\n", "\n", "One very important point of consideration is that xarray data variables (in both Datasets and DataArrays) can store both unit-aware and unit-naive array types. Unit-naive array types will be used by default in xarray, so we need to convert to a unit-aware type if we want to use xarray operations while preserving unit correctness. MetPy provides the .quantify() method for this (named since we are turning the data stored inside the xarray object into a Pint Quantity object)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'Temperature_isobaric' (time1: 9, isobaric3: 31,\n",
       "                                          latitude: 81, longitude: 131)>\n",
       "<Quantity([[[[258.      258.      257.9     ... 258.2     258.2     258.3    ]\n",
       "   [257.9     257.9     257.8     ... 258.2     258.2     258.2    ]\n",
       "   [257.8     257.8     257.8     ... 258.1     258.1     258.1    ]\n",
       "   ...\n",
       "   [260.8     260.8     260.8     ... 258.      258.1     258.2    ]\n",
       "   [261.      261.      261.      ... 258.1     258.1     258.2    ]\n",
       "   [261.2     261.2     261.2     ... 258.1     258.2     258.3    ]]\n",
       "\n",
       "  [[256.5     256.5     256.5     ... 255.6     255.6     255.5    ]\n",
       "   [256.4     256.4     256.3     ... 255.5     255.5     255.5    ]\n",
       "   [256.2     256.2     256.2     ... 255.4     255.4     255.4    ]\n",
       "   ...\n",
       "   [259.1     259.1     259.1     ... 261.1     261.4     261.6    ]\n",
       "   [259.2     259.2     259.2     ... 261.1     261.4     261.6    ]\n",
       "   [259.4     259.3     259.3     ... 261.1     261.4     261.6    ]]\n",
       "\n",
       "  [[251.4     251.4     251.4     ... 250.5     250.5     250.5    ]\n",
       "   [251.2     251.3     251.3     ... 250.4     250.4     250.4    ]\n",
       "   [251.1     251.1     251.1     ... 250.3     250.3     250.2    ]\n",
       "   ...\n",
       "...\n",
       "   ...\n",
       "   [295.40698 295.40698 295.40698 ... 295.70697 295.70697 295.507  ]\n",
       "   [295.30698 295.30698 295.30698 ... 295.20697 295.30698 295.40698]\n",
       "   [295.30698 295.107   295.107   ... 294.80698 295.007   295.107  ]]\n",
       "\n",
       "  [[291.9     290.2     290.5     ... 288.4     288.6     288.9    ]\n",
       "   [291.4     290.7     290.8     ... 289.3     289.5     289.9    ]\n",
       "   [290.7     290.4     290.5     ... 289.9     290.1     290.5    ]\n",
       "   ...\n",
       "   [297.5     297.6     297.5     ... 297.9     297.8     297.7    ]\n",
       "   [297.4     297.4     297.3     ... 297.4     297.5     297.6    ]\n",
       "   [297.2     297.1     297.      ... 297.      297.2     297.3    ]]\n",
       "\n",
       "  [[293.3     291.6     291.9     ... 286.8     286.9     287.     ]\n",
       "   [292.8     292.1     292.2     ... 287.      286.9     287.     ]\n",
       "   [292.1     291.8     291.9     ... 287.1     287.3     287.5    ]\n",
       "   ...\n",
       "   [299.6     299.7     299.7     ... 300.1     300.      299.8    ]\n",
       "   [299.5     299.5     299.4     ... 299.5     299.7     299.7    ]\n",
       "   [299.4     299.3     299.2     ... 299.2     299.3     299.4    ]]]], 'kelvin')>\n",
       "Coordinates:\n",
       "  * time1      (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...\n",
       "    reftime    datetime64[ns] 2017-09-05T12:00:00\n",
       "  * latitude   (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n",
       "  * longitude  (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n",
       "    metpy_crs  object Projection: latitude_longitude\n",
       "  * isobaric3  (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05\n",
       "Attributes:\n",
       "    long_name:                      Temperature @ Isobaric surface\n",
       "    Grib_Variable_Id:               VAR_0-0-0_L100\n",
       "    Grib2_Parameter:                [0 0 0]\n",
       "    Grib2_Parameter_Discipline:     Meteorological products\n",
       "    Grib2_Parameter_Category:       Temperature\n",
       "    Grib2_Parameter_Name:           Temperature\n",
       "    Grib2_Level_Type:               100\n",
       "    Grib2_Level_Desc:               Isobaric surface\n",
       "    Grib2_Generating_Process_Type:  Forecast\n",
       "    grid_mapping:                   LatLon_361X720-0p25S-180p00E
" ], "text/plain": [ "\n", "\n", "Coordinates:\n", " * time1 (time1) datetime64[ns] 2017-09-05T12:00:00 ... 2017-09-06T12:0...\n", " reftime datetime64[ns] 2017-09-05T12:00:00\n", " * latitude (latitude) float32 50.0 49.5 49.0 48.5 ... 11.5 11.0 10.5 10.0\n", " * longitude (longitude) float32 250.0 250.5 251.0 251.5 ... 314.0 314.5 315.0\n", " metpy_crs object Projection: latitude_longitude\n", " * isobaric3 (isobaric3) float64 100.0 200.0 300.0 ... 9.5e+04 9.75e+04 1e+05\n", "Attributes:\n", " long_name: Temperature @ Isobaric surface\n", " Grib_Variable_Id: VAR_0-0-0_L100\n", " Grib2_Parameter: [0 0 0]\n", " Grib2_Parameter_Discipline: Meteorological products\n", " Grib2_Parameter_Category: Temperature\n", " Grib2_Parameter_Name: Temperature\n", " Grib2_Level_Type: 100\n", " Grib2_Level_Desc: Isobaric surface\n", " Grib2_Generating_Process_Type: Forecast\n", " grid_mapping: LatLon_361X720-0p25S-180p00E" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "temperature = data_parsed['Temperature_isobaric'].metpy.quantify()\n", "temperature" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how the units are now represented in the data itself, rather than as a text attribute. Now, even if we perform some kind of xarray operation (such as taking the zonal mean), the units are preserved." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "tags": [] }, "outputs": [], "source": [ "# Take the mean over time at the 1000 hPa level\n", "temperature.sel(isobaric3=1000).mean('time1')\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Plotting \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Load data\n", "ds = xr.open_dataset(get_test_data('narr_example.nc', as_file_obj = False))\n", "\n", "# Parse full dataset\n", "ds = ds.metpy.parse_cf()\n", "\n", "# Grab lat/lon values from file as unit arrays\n", "lats = ds.lat.metpy.unit_array\n", "lons = ds.lon.metpy.unit_array\n", "\n", "# Get the valid time\n", "vtime = ds.Temperature.metpy.time[0]\n", "\n", "# Get the 700-hPa heights without manually identifying the vertical coordinate\n", "hght_700 = ds.Geopotential_height.metpy.sel(vertical=700 * units.hPa,\n", " time=vtime)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Cartopy CRS, useful for plotting\n", "hght_700.metpy.cartopy_crs" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import cartopy for plotting\n", "import cartopy.crs as ccrs\n", "import cartopy.feature as cfeature\n", "import matplotlib.pyplot as plt\n", "\n", "# Open figure\n", "fig = plt.figure(figsize=(12, 12))\n", "\n", "# Add a plot with coordinate axes (crs)\n", "ax = fig.add_subplot(1, 1, 1, projection=hght_700.metpy.cartopy_crs)\n", "\n", "# Get coordinates\n", "x = ds.x\n", "y = ds.y\n", "\n", "# Plot data\n", "ax.imshow(hght_700, extent=(x.min(), x.max(), y.min(), y.max()),\n", " cmap='RdBu', origin='lower' if y[0] < y[-1] else 'upper')\n", "\n", "# Add coastlines\n", "ax.coastlines(color='tab:blue', resolution='10m')\n", "\n", "# Show the plot\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculations" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from metpy.units import units\n", "import metpy.calc as mpcalc" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "temperature = [20] * units.degC\n", "rel_humidity = [50] * units.percent\n", "print(mpcalc.dewpoint_from_relative_humidity(temperature, rel_humidity))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "speed = np.array([5, 10, 15, 20]) * units.knots\n", "direction = np.array([0, 90, 180, 270]) * units.degrees\n", "u, v = mpcalc.wind_components(speed, direction)\n", "print(u, v)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Reminder \n", "
    \n", "
  • Import the package, aka import metpy\n", "
  • Metpy works with units\n", "
  • MetPy aims to have three primary purposes: read and write meteorological data (I/O), calculate meteorological quantities with well-documented equations, and create publication-quality plots of meteorological data\n", "
\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "Homework: Check out the metpy tutorial and get familiar with the syntax.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "meteo", "language": "python", "name": "meteo" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.5" } }, "nbformat": 4, "nbformat_minor": 4 }