<jupyter_start><jupyter_text>Introduction to SatpySatpy is another third-party Python library that provides additional functionality for specifically working with satellite instrument data. Satpy tries to ease the common tasks that people have to perform when working with satellite data. This lesson will provide an introduction to Satpy for reading satellite imagery data from the GOES-16 ABI instrument, creating color RGB images, and plotting data using matplotlib.This lesson is a shortened version of the lessons available in the Pytroll Satpy Half-day tutorial. For more details on Satpy or a more thorough explanation of the topics discussed in this lesson, view that tutorial [here](https://github.com/pytroll/tutorial-satpy-half-day). What are we working with?We will be working with earth-observing meteorological satellite instrumentdata. There are many different variations of how the data can be structured,what it physically represents, how it changes over time, and how it can be usedin a particular type of application. The following is an overview of some commoncharacteristics of meteorological satellite observations. We'll go in to a fewdetails as we explore real data later on and completely ignore otherdetails for simplicity. Earth-observingSatpy operates on data from earth-observing satellite instruments. The datacan be used to study changes in the atmosphere, vegetation, oceans,pollution, and many others. Geostationary versus Polar-orbitingSatellites can be in a high-altitude geostationary orbit or a lower altitudepolar-orbiting orbit. Geostationary satellites typically provide observationsat higher temporal rates of the same region faster (~30s in some), but usuallyhave a lower spatial resolution and only see part of the Earth.Polar orbiters are usually at a higher spatial resolution and covermuch more of the Earth, but it takes longer to get all of this coverage.            Geostationary        Polar-orbiting                                                                Credit: <NAME>, SSEC, UW-Madison ImageryThe instrument data we will be working with is imagery data; data from"imager" instruments. In most cases, these arrays of data points can bethought of as a 2D image of pixels.                                        &rarr;                                 GeolocatedThe data we will be looking at is geolocated. We need to be able to assign eachpixel of data to a geographic region. Data footprints that can be somewhatdifficult to describe are typically simplified by specifyingonly the center point and the pixel's radius or cell width. We may know theexact longitude and latitude coordinates of each pixel or we may be given agridded version of the data where each pixel is spread uniformly across arectangular area.If you are familiar with projections, we'll get to those later. BandsThe satellite instrument data we will be working with are a collectionof radiation measurements. By taking individual measurements ofthe radiation reflected or emitted by objects on Earth, we can get a goodsnapshot of the Earth from space. Satellite instruments will typicallyhave multiple bands or channels where each one measures a specific wavelengthof the electromagnetic spectrum. Each band can show us something differentabout the Earth. SectorsSome geostationary satellites record data for different sectors or regionsof the Earth depending on their purpose and design. For example the GOES-16ABI instrument has 4 sectors:- Full Disk (FD or FLDK)- Continental United States (CONUS)- Mesoscale 1 (M1)- Mesoscale 2 (M2)The below image shows these different sectors and how quickly GOES-16 ABI records data for them. The colors on the image indicate how much of the Earth each data pixel measures.Credit: <NAME>, CIMSS Reading data with the `Scene`Satpy's main interface for working with data is the `Scene` class. We can provide the `Scene` with data files and load them with a "reader". In this notebook we'll explore the basic data loading and data access functionality provided by Satpy while also providing a basic introduction to xarray's `DataArray` objects and `dask` arrays.Before importing and using Satpy, we run some python code to do some initial setup. This includes turning off warnings and limiting the number of resources we use. These are precautions to make these examples work on the most machines.<jupyter_code>%run ../init_notebook.py
from satpy import Scene
from glob import glob

# Get the list of GOES-16 ABI files to open
filenames = glob('../data/abi_l1b/20180511_texas_fire_abi_l1b_conus/*.nc')
len(filenames)
scn = Scene(reader='abi_l1b', filenames=filenames)<jupyter_output><empty_output><jupyter_text>We've now created a `Scene` object. Under the hood Satpy has sorted the files and determined what we can access. We haven't actually loaded any data so our dict-like `Scene` object is empty. To find out what data can be loaded from the file we can use the `available_dataset_names`.<jupyter_code>scn.available_dataset_names()<jupyter_output><empty_output><jupyter_text>The `Scene` is telling us that we have all 16 ABI channels available to load. This list includes any product that we can load from the file that the "abi_l1b" reader is configured to access. If we didn't provide all of the necessary files or the data was missing from the file for some reason, that product would not be listed here.| Channel     | Wavelength  |  Resolution  || ----------- | ----------- |  ----------- || C01         | 0.47µm      |  1000m       || C02         | 0.64µm      |  250m        || C03         | 0.64µm      |  1000m       || C04         | 1.37µm      |  2000m       || C05         | 1.60µm      |  1000m       || C06         | 2.20µm      |  2000m       || C07         | 3.90µm      |  2000m       || C08         | 6.20µm      |  2000m       || C09         | 6.90µm      |  2000m       || C10         | 7.30µm      |  2000m       || C11         | 8.40µm      |  2000m       || C12         | 9.60µm      |  2000m       || C13         | 10.30µm     |  2000m       || C14         | 11.20µm     |  2000m       || C15         | 12.30µm     |  2000m       || C16         | 13.30µm     |  2000m       |Let's pick one of these channels, load it, and look what information is provided by Satpy.<jupyter_code>my_channel = 'EDITME'
scn.load([my_channel])
# use brackets to access products like a normal dict
scn[my_channel]<jupyter_output><empty_output><jupyter_text>Xarray and DaskAbove we see an `xarray.DataArray` object with a lot of metadata.There are a few elements to get familiar with when working with DataArray's from Satpy:* `dask.array`: We don't see any actual imagery data. Our data is stored in a `dask` array instead of a traditional numpy array. This means our data's loading and calculations are delayed.* `Attributes`: A dictionary where the metadata is stored. Some is from the file, some is added by the "abi_l1b" reader to assist future Satpy operations. Some of the more important keys are:  * `platform_name`  * `sensor`  * `name`  * `wavelength`  * `units`  * `calibration`  * `standard_name`  * `start_time`  * `area` (more on this later)If we want to access the attributes, we use the `.attrs` attribute.<jupyter_code>scn[my_channel].attrs['start_time']<jupyter_output><empty_output><jupyter_text>We can access the dimension names of the data using the `.dims` attribute.<jupyter_code>scn[my_channel].dims<jupyter_output><empty_output><jupyter_text>The sizes of those dimensions:<jupyter_code>scn[my_channel].sizes['y']<jupyter_output><empty_output><jupyter_text>DataArrays also provide us access to traditional numpy properties like `shape` and `ndim`.<jupyter_code>scn[my_channel].shape
scn[my_channel].ndim<jupyter_output><empty_output><jupyter_text>Use them like numpy arraysIn most cases, Xarray's `DataArray` objects can be used just like a regular numpy array. When the actual data values are needed they will be computed. This allows us to use `DataArray` objects with other python tools with little to no extra work. Then do a simple matplotlib plot to view our data.Note: If running on a Jupyter Lab session you may need to change "notebook" in the below cell to "inline".<jupyter_code>%matplotlib widget

import matplotlib.pyplot as plt
plt.figure()
plt.imshow(scn[my_channel])
plt.colorbar()<jupyter_output><empty_output><jupyter_text>We can use matplotlib calls manually like above, but Xarray also provides its own plotting utility functions to make this easier.<jupyter_code>plt.figure()
scn[my_channel].plot.imshow(cmap='viridis')<jupyter_output><empty_output><jupyter_text>If we compare the metadata in our DataArray's `.attrs` with what labels are on the plot, we can see where Xarray has made its best guess about what to name components of the plot. It used attributes like `long_name` for the colorbar and the names of the dimensions for the axis labels. Xarray's plotting utilities are simple wrappers around matplotlib so we still have access to everything from matplotlib. We can add common matplotlib function calls like `plt.title(my_channel)` to the above cell, for example, to change the title.We can also change the colormap by passing the `cmap` keyword argument to the call to `imshow` (ex. `cmap='viridis'`). For a full list of the builtin matplotlib colormaps see the [matplotlib documentation](https://matplotlib.org/tutorials/colors/colormaps.html). By default matplotlib will use `viridis` but we can also try others like `plasma`, `magma`, `RdBu_r`, `Reds`, or `tab20b`.  Exercise**Time: 5-10 minutes**Using the above examples as a guide, load additional channels, view them with matplotlib, and explore the data by using matplotlib's interactive notebook widget (see toolbar at the bottom left of the image).<jupyter_code># Add your code here<jupyter_output><empty_output><jupyter_text>RGBs and Other CompositesMeteorological satellite instrument bands provide a look at the Earth at various wavelengths. Each band represents the electro-magnetic radiation observed at one wavelength or range of wavelengths. While this is very useful, flexible, and requires a very limited amount of processing to visualize, it does require that the viewer know what to expect from each band or wavelength. In addition to viewing individual bands, another common option for viewing satellite data is to combine multiple bands in to one image such as RGB (Red, Green, Blue) images. By assigning certain bands to particular channels in the output RGB image we can provide an easy way to identify features in the data. Instead of having to know what a particular wavelength is sensitive to, we can now focus on what colors in particular RGB recipe mean.Using our ABI CONUS data from before, let's look at some of the built-in RGB composites provided by Satpy.<jupyter_code>filenames = glob('../data/abi_l1b/20180511_texas_fire_abi_l1b_conus/*.nc')
scn = Scene(reader='abi_l1b', filenames=filenames)
scn.available_composite_names()<jupyter_output><empty_output><jupyter_text>The above list shows what composites Satpy knows how to make and that it has the inputs for. A lot of the recipes in Satpy come from standard recipes from science agencies like EUMETSAT, NOAA, or NASA. Let's first look at the `airmass` RGB. The `airmass` RGB is made of the following bands:```R: C08 - C10G: C12 - C13B: C08```The red channel is the difference between the C08 (6.185µm) and C10 (7.34µm) bands, the green channel is the difference between the C12 (9.61µm) and C13 (10.35µm) bands, and the blue channel is the C08 (6.185µm) band.The EUMETSAT EUMeTrain manual summarizes the use of the `airmass` RGB as:> As both water vapor channels (WV 6.2 µm and WV 7.3 µm) are included in this composite, the main applications are the detection of dynamic processes, such as rapid cyclogenesis, jet streams and PV anomalies.Read more on the [EUMeTrain website](http://www.eumetrain.org/data/3/306/navmenu.php?tab=9&page=4.0.0).In Satpy, we can load preconfigured composites like any normal file-provided product. Satpy will look at the recipe, determine what bands need to loaded, what modifications/corrections need to be applied to satisfy the recipe, and then give us the `DataArray` representing that RGB back.<jupyter_code>scn.load(['airmass'])
scn['airmass']<jupyter_output><empty_output><jupyter_text>The key part of the above xarray output is that this is the first 3 dimensional DataArray we've worked with in Satpy with the `(bands: 3, y: 1500, x: 2500)` dimensions. The `bands` dimension tells us what each "band" represents. In this case Red, Green, and Blue channels of an RGB image. If we want to display our data in a matplotlib figure we need to convert the values in to something matplotlib understands. Right now we have the real observed values from the satellite instrument, we need to scale them to fit in a traditional image format like PNG. Satpy provides a `get_enhanced_image` utility function to simplify this.We can tell xarray what dimension of our `DataArray` represents the "RGB" axis. We do this by passing the `rgb` keyword argument with the name of the dimension. In Satpy this dimension should always be called `'bands'`. Without this matplotlib would not be able to plot our RGB properly.Lastly, we'll tell the `imshow` method that we *know* that our RGB will be between 0 and 1 with the `vmin` and `vmax` keyword arguments so there is no need to compute the min/max of the data.Let's plot our `airmass` RGB composite.<jupyter_code>%matplotlib widget
import matplotlib.pyplot as plt
from satpy.writers import get_enhanced_image

fig = plt.figure()
img = get_enhanced_image(scn['airmass'])
# get DataArray out of `XRImage` object
img_data = img.data
img_data.plot.imshow(rgb='bands', vmin=0, vmax=1)<jupyter_output><empty_output><jupyter_text>Compositing Different ResolutionsWhen generating the `airmass` RGB we were fortunate in that all of the input data was the same resolution. Let's look at the `natural_color` RGB now, where the input resolutions differ. The `natural_color` RGB is useful as described by EUMeTrain for seeing snow and ice:> Snow on the ground as well as ice over mountains, frozen lakes and sea ice appear cyan in the Natural Colour RGB images.See the [EUMeTrain website](http://www.eumetrain.org/resources/MSG_natcolour_cyan_snow.html) for more information.For the ABI instrument this composite is a combination of the following channels:```R: C05 (1000m)G: C03 (1000m)B: C02 (500m)```This means that we are trying to combine arrays of different shapes. The easiest way to make data the same resolution and make it easier to compare or combine, is to resample. Let's load this composite and see how Satpy handles generating a composite like this.<jupyter_code>scn.load(['natural_color'])
scn['natural_color']<jupyter_output><empty_output><jupyter_text>We didn't do anything wrong, but Satpy doesn't have a `natural_color` composite in the `Scene`. Satpy knows that the reason it wasn't able to generate this composite was because the data was not all the same resolution. The `Scene` knows that we need to resample before we can be able to generate this composite. Let's do that now by using the `native` resampler. We'll also use the `min_area` method to aggregate the data to the lowest resolution in the Scene; 2km in this case (from the `airmass` dependencies):<jupyter_code>new_scn = scn.resample(scn.min_area(), resampler='native')
new_scn['natural_color']<jupyter_output><empty_output><jupyter_text>Once resampling had completed, the `Scene` attempted to regenerate any composites that were previous requested. Once successfully generated the `natural_color` RGB shows up like any other product. Let's follow the same procedure we used for plotting the `airmass` RGB, but for the `natural_color` RGB. Note we are using the `new_scn` because the `natural_color` RGB only exists in this resampled `Scene`.<jupyter_code>from dask.diagnostics import ProgressBar
with ProgressBar():
    plt.figure()
    img = get_enhanced_image(new_scn['natural_color'])
    # get DataArray out of `XRImage` object
    img_data = img.data
    img_data.plot.imshow(vmin=0, vmax=1, rgb='bands')<jupyter_output><empty_output><jupyter_text>Saving imagesLooking at satellite instrument data in a Notebook is nice, but it is difficult for Notebook widgets to provide scientists with everything they need for analyzing a dataset. There are many applications that exist for viewing images on a map that provide a suite of tools for analyzing and manipulating the data. One way to get data (images) in to these applications is by providing it in a special image format called GeoTIFF. Luckily, Satpy makes it very easy to save any Scene to a series of GeoTIFFs by using the `save_datasets` method. Let's recreate out `Scene` object and call this method.<jupyter_code>scn = Scene(reader='abi_l1b', filenames=filenames)
scn.load(['C07'])
with ProgressBar():
    scn.save_datasets()<jupyter_output><empty_output><jupyter_text>Unfortunately, if you are running this tutorial from a remote BinderHub then you won't be able to display the GeoTIFF since it isn't a supported image format by most browsers. You can change the format that Satpy creates. Instead of a GeoTIFF (`.tif`) we can make a PNG (`.png`) by specifying the writer parameter:<jupyter_code>with ProgressBar():
    scn.save_datasets(writer='simple_image')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.06392168154333429, 'fraction_numerical': 0.01618197523754679, 'mean_word_length': 4.449011609664261, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 77, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '116399', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4965, 'n_tokens_neox': 4579, 'n_words': 2296, 'path': 'notebooks/02_intro_to_satpy.ipynb', 'repo_name': 'djhoese/tutorial-intro-python-satpy', 'size': '17381'}