<jupyter_start><jupyter_text>Satpy - GOES-16 ABI Example**Author:** <NAME>In this notebook, we use the python [Satpy](https://satpy.readthedocs.io/en/latest/) library to read from NOAA's GOES-16 ABI instrument, resample data, and write it back out to a file. The Satpy library takes advantage of many open source libraries:- xarray and dask for data containers and optimization- NetCDF4, HDF5, and rasterio libraries for file reading and writing- Pyresample library for resamplingWe get our data from the [publicly available data](https://console.cloud.google.com/storage/browser/gcp-public-data-goes-16) NOAA has put on Google Cloud Platform (GCP). We'll be specifically using data from March 13th, 2019 (March 14th UTC) over the United States to look at an intense midlatitude cyclone. This case was chosen from a blog post by <NAME> on the CIMSS Satellite Blog. Scott's blog post can be found [here](http://cimss.ssec.wisc.edu/goes/blog/archives/32345) which includes much more in-depth analysis.We'll start by importing the various python tools we'll use and some basic configuration to limit how much of our machine we use.<jupyter_code>import os
from glob import glob
import urllib.request
import warnings
warnings.filterwarnings('ignore')
os.environ['OMP_NUM_THREADS'] = "2"
os.environ['PYTROLL_CHUNK_SIZE'] = "1024"

import dask
from dask.diagnostics import ProgressBar
dask.config.set(num_workers=4)
import gcsfs
from satpy import Scene, demo<jupyter_output><empty_output><jupyter_text>Download example dataFirst we download the GOES-16 ABI NetCDF files we will be working with. There are 16 bands on the ABI instrument so we download 16 files for a single time step of data. The data is over the "CONUS", or Continental United States, region. We use Satpy's `demo` module, using the `gcsfs` package underneath, to download these files from the public GCP bucket storage.<jupyter_code>filenames = demo.get_us_midlatitude_cyclone_abi()<jupyter_output><empty_output><jupyter_text>Load data with SatPyNow that we have put our data somewhere we can access it, we pass the files to SatPy to analyze and open. We create a `Scene` object which acts as a container for our data as we load and analyze it. After creating the `Scene` we ask what channels (a.k.a. bands) we can load from the files.<jupyter_code>scn = Scene(reader='abi_l1b', filenames=filenames)
scn.available_dataset_names()<jupyter_output><empty_output><jupyter_text>We can also create composite images by combining one or more bands together. Here we ask, based on the available channels, what composites can be make:<jupyter_code>scn.available_composite_names()<jupyter_output><empty_output><jupyter_text>To actually load these bands or create these composites we use the `load` method. Here we ask for all 16 channels as well as a `true_color`, `natural_color`, and `airmass` RGB:<jupyter_code>all_channels = ['C{:02d}'.format(x) for x in range(1, 17)]
composites = ['true_color', 'natural_color', 'airmass']
scn.load(all_channels + composites)<jupyter_output><empty_output><jupyter_text>The `Scene` object now acts as a dictionary-like container of our data. Each band or "product" that we've loaded can be accessed by it's name **or** by its wavelength in the case of the bands. The value that gets returned is an `xarray.DataArray` object.<jupyter_code>scn['C01']
# scn[0.47]
scn['airmass']<jupyter_output><empty_output><jupyter_text>Although some of the RGBs, like "airmass", were created completely from the `load` call, some could not be generated yet. Usually this means that further processing is needed to get the various RGB inputs in to a compatible state. The most common case of this is bands being at different resolutions; typically resolved by resampling data. Let's investigate the current state of the Scene where we see that "true_color" and "natural_color" are missing:<jupyter_code>scn.missing_datasets<jupyter_output><empty_output><jupyter_text>Resample dataA common operation done with satellite data is to resample the data to a gridded area with uniform pixel sizes. In addition to helping us resolve bands to the same resolution (see above), this also allows us to change the geographic projection of the data or to crop out a specific region of interest to study. SatPy currently provides a "native" resampling method for the cases where we don't want to change the projection or coverage, but only unify the resolutions of the currently loaded data. This is fast because we are only replicating or aggregating data points. There is also a "nearest" or nearest neighbor resampling method to handle the other cases.We start by using the "native" resampler to make all bands match the resolution of the lowest resolution band (2km) to save memory and processing time. This call will automatically detect that "true_color" and "natural_color" were not created, but can be now. We'll use this Scene later on to make full resolution cloud-optimized geotiffs.<jupyter_code>new_scn = scn.resample(scn.min_area(), resampler='native')
new_scn.missing_datasets<jupyter_output><empty_output><jupyter_text>Saving to diskNow we want to save our data to cloud-optimized geotiffs (COGs). These are normal geotiffs that are structured in a way to make them easier to read from cloud-computing sources and environments. See [COG](https://www.cogeo.org/) for more information.We add the optional `ProgressBar` context manager from dask to get an idea of where we are in the processing. Due to the way SatPy takes advantage of dask, we haven't really done any number crunching until this point. The `save_datasets` method below will write these geotiffs to the current directory or we could use the `base_dir` option to specify an output directory. We'll use the `datasets` keyword argument to limit what data we write out as geotiffs.<jupyter_code>with ProgressBar():
    new_scn.save_datasets(datasets=['C01', 'airmass'], compress='DEFLATE', tiled=True, copy_src_overviews=True)<jupyter_output><empty_output><jupyter_text>Cropping dataSimilar to resampling to a different projection or a different region, we can also perform the simpler operation of cropping the data to a specific bounding box. This will slice the data to a smaller size and make it faster to perform more complex analysis on a specific region we're curious about. Since our new smaller region is relatively small to our original data, we can use the `show` method without worrying too much about making an image too large for our browser.<jupyter_code>cropped_scn = new_scn.crop(ll_bbox=(-103.0, 35.0, -99.0, 41.0))
print("Uncropped: ", new_scn['true_color'].shape)
print("Cropped: ", cropped_scn['true_color'].shape)<jupyter_output><empty_output><jupyter_text>Now we'll use Xarray's plotting utilities to create a matplotlib figure of our data. Right now the 'true_color' image still has data values corresponding to the input bands (i.e. reflectance percentages). We will use Satpy's `get_enhanced_image` function to scale our data values to something that can be displayed as an image. This takes about 20-30 seconds to generate due to all of the atmospheric corrections being performed to generate the image.<jupyter_code>%matplotlib inline

from satpy.writers import get_enhanced_image
# normalize the data to be shown as an image
enhanced_img = get_enhanced_image(cropped_scn['true_color']).data

with ProgressBar():
    # matplotlib needs (rows, columns, 3)
    # persist data in memory so it is only computed once
    enhanced_img = enhanced_img.transpose('y', 'x', 'bands').persist()
# use xarray's plotting utilities
enhanced_img.plot.imshow()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07022140704834272, 'fraction_numerical': 0.00877767588104284, 'mean_word_length': 5.644038294168842, 'pattern_counts': {'":': 0, '<': 47, '<?xml version=': 0, '>': 47, 'https://': 3, 'lorem ipsum': 0, 'www.': 1, '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': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '296983', 'license': 'CC0-1.0', 'n_tokens_mistral': 2229, 'n_tokens_neox': 2122, 'n_words': 983, 'path': 'satpy_earth_observations_goes16_abi.ipynb', 'repo_name': 'pangeo-data/open-source-geoscience', 'size': '7647'}