Subsetting NSIDC Data in the Cloud and Exporting to Your Local Machine

1. Tutorial Overview

This notebook demonstrates searching for and directly accessing cloud-hosted granules from an Amazon Compute Cloud (EC2) instance using the earthaccess package. It walks through subsetting the data by variable, clipping it spatially to an area of interest using a provided outline, and finally copying the processed files to your local machine.

Direct Access is an efficient way to work with data stored in an S3 Bucket when you’re working in the cloud. Cloud-hosted granules can be opened and loaded into memory without the need to download them first. This means you can subset and clip data before it ever touches your local storage, keeping the volume hitting your local machine to a minimum.

Many NASA earthdata collections offer customization services through Harmony, which allows for subsetting, reprojecting, and/or reformatting data before you download or open it. Please see the following tables for ICESat-2 and SMAP for which NSIDC datasets provide these customization services.

However, as seen from the tables linked above, there are many data sets that do not yet have customization services. For these datasets, it is then necessary to either download the full sized granules, or to develop a workflow to manually customize granules in the cloud before transferring them to local storage. This tutorial was developed to walk through the latter option.

For this tutorial, we’ll use granules from the High Mountain Asia 12 km Modeled Estimates of Aerosol Transport, Chemistry, and Deposition Reanalysis, 2003-2019, Version 1 (HMA2_MATCHA) data set. HMA2_MATCHA contains a 12 km resolution, simulated reanalysis of aerosol transport, chemistry, and deposition over the High Mountain Asia (HMA) region for 1 January 2003 through 31 August 2019. Also known as the Model for Atmospheric Transport and Chemistry in Asia (MATCHA), the data comprise a wide range of variables intended to help assess the impacts of aerosols on the cryosphere in the HMA region. Data are stored in NetCDF-4 format.

Learning Objectives

By the end of this demonstration you will be able to:
1. use earthaccess to search for and directly access cloud-hosted HMA2_MATCHA data in an AWS EC2 instance
2. subset the dataset to specific variables using xarray
3. spatially clip the dataset to an area of interest using rioxarray
4. export the variable and spatially subsetted dataset and transfer it to your local machine

Prerequisites

  1. An AWS account with an EC2 instance in the United States (Oregon) — us-west-2 region.
    • A step-by-step article for setting up an AWS EC2 instance and installing Conda is provided here.
    • This tutorial was developed and tested on an AWS EC2 t2.medium Ubuntu instance with a 20 GB storage volume.
  2. An Earthdata Login account is required for data access. If you don’t have one, you can register for one here.
  3. A .netrc file, that contains your Earthdata Login credentials, in the home directory of your EC2 instance. Instructions for creating a .netrc are provided here.
  4. This notebook file on your EC2 instance.
    • This notebook can be downloaded directly from the NSIDC Data Cookbook Repository to your EC2 instance using the command:
      REPLACE WITH MAIN BRANCH
    wget https://raw.githubusercontent.com/nsidc/NSIDC-Data-Cookbook/refs/heads/subsetting-tutorial/tutorials/HMA_subsetting_cloud_workflow.ipynb
  5. An environment with all the necessary modules installed on your EC2 instance.
    • An environment.yml file is available in the root of the NSIDC Data Cookbook repository. This file can be used to create the conda environment on your EC2 instance with the command:
    conda env create -f https://raw.githubusercontent.com/nsidc/NSIDC-Data-Cookbook/refs/heads/main/environment.yml

2. Set up

Import Packages

# For searching NASA data
import earthaccess

# For reading, processing, and plotting data
import xarray as xr
import rioxarray as rio
import hvplot.xarray
import geopandas as gpd
from affine import Affine

# For crafting output filenames
from pathlib import Path

# For memory management
import gc

Authenticate

The first step is to get the correct authentication that will allow us to get cloud-hosted HMA data. This is all done through Earthdata Login. The login method also gets the correct AWS credentials.

Login requires your Earthdata Login username and password. The login method will automatically search for these credentials as environment variables or in a .netrc file, and if those aren’t available it will prompt us to enter our username and password. We use a .netrc strategy. A .netrc file is a text file located in our home directory that contains login information for remote machines. If we don’t have a .netrc file, login can create one for us.

auth = earthaccess.login()

Search for granules using earthaccess.search_data() with temporal and granule filters

The HMA2_MATCHA dataset provides two-dimensional surface (SFC) data at one hour intervals, and three-dimensional atmospheric (ATM) data at three-hour intervals. For this tutorial, we’ll filter the results to only return atmospheric data granules. To do this, we use the granule_name argument with a wildcard search (*) to only return granule names containing “*_ATM_*”.

We will print the number of results returned using the chosen filters.
For fun, we’ll also print the total storage that would be needed if we downloaded the results as-is to our local machine.

results = earthaccess.search_data(
    short_name='HMA2_MATCHA',
    cloud_hosted=True,
    temporal=('2018-07-01', '2018-07-01'),
    granule_name='*_ATM_*',
)

print(f"Number of results: {len(results)}")

results_volume = sum(result.size() for result in results)
print(f"Total volume of results: {round(results_volume, 1)} MBs ({round(results_volume / 1000, 1)} GBs)")

Use direct access to open the data stored on S3

The files can now be opened using the open method, with the auth object created at the start of the notebook providing the Earthdata Login authentication and AWS credentials.

files = earthaccess.open(results)

3. Subset the dataset to specific variables

Open the file and explore the variables

We will start with viewing and processing just a single file from the results. To open the file, we use decode_coords='all' to load the grid mapping variable (crs, in this case) into the data set coordinates, and mask_and_scale=False to ensure all original data values are maintained.

The structure of the NetCDF files for HMA2_MATCHA are relatively simple, with all the variables being located in the root of the file. Because there are no variables contained within nested groups, we can use xarray’s open_dataset() method to open the file and have access to all the variables in the file.

ds = xr.open_dataset(files[0], engine='h5netcdf', decode_coords='all', mask_and_scale=False)
ds

Create a variables list and subset the dataset

We can see the HMA2_MATCHA atmospheric files contain 128 data variables. The attributes for each variable in the file can be explored by clicking the ‘Data Variables’ dropdown in the above output. This can be useful for initial exploration of the dataset and for determining which variables may be useful to our specific analysis.

For this tutorial we want the data variables lat (latitude), lon (longitude), time, temperature, ua (x wind component), va (y wind component), and bc_3d_tot (total black carbon mass concentration), which we store as strings to the data_vars list below.

  • Note: The x, y, and lev dimensions are used by our chosen data variables, so these dimensions will transfer to our subsetted dataset automatically and don’t need to be specified.
data_vars = [
    'lat',
    'lon',
    'time',
    'temperature',
    'ua',
    'va',
    'bc_3d_tot'
]

With the data_vars list defined, we can create a subset of the original dataset that contains only those variables:

ds_subset = ds[data_vars]
ds_subset
  • Note: We included time in the data_vars list, but the output above only shows time as a dimension and coordinate variable – not a data variable. In the original dataset, time is a coordinate variable, but it is not associated with any of the data variables in our data_vars list. As a result, it would be dropped during subsetting unless we explicitly include it. Because the time coordinate contains important information, we include it in data_vars to ensure it is retained.

Let’s extract the uncompressed size of the original dataset, the subsetted data, and compare the savings on volume.

ds_volume = ds.nbytes
ds_subset_volume = ds_subset.nbytes
print(f"The beginning volume of the file (uncompressed) was {round((ds_volume / 1e+6), 1)} MBs")
print(f"The volume of the subsetted file (uncompressed) is now {round((ds_subset_volume / 1e+6), 1)} MBs")
print(f"That is a {round(100 - ((ds_subset_volume / ds_volume)*100), 1)} percent reduction in volume for this one file.")

Plot a variable of the subsetted dataset

We can now use hvplot to plot the bc_3d_tot variable of our subsetted data.

def plot_variable(ds, var):
    return ds[var].hvplot.quadmesh(
        x='x',
        y='y',
        groupby='lev',
        project=True,
        tiles=True,
        geo=True,
        rasterize=True,
        cmap='turbo',
)

plot_variable(ds_subset, 'bc_3d_tot')

4. Clip the dataset to your area of interest using a provided outline

Load the outline for your area of interest

From the plot above we can see the HMA2_MATCHA dataset covers the whole High-Mountain Asia region. For this tutorial, we want to clip this down to only include data over the Himalayan mountain range.

Let’s use GeoPandas (aliased as gpd) to load and plot the outline of the Himalayas. We’ll use the Himalayas_outline.geojson file that lives in the example_data/subsetting_tutorial directory of the NSIDC Data Cookbook repository.
We can use the GitHub raw URL with GeoPandas to load the GeoJSON file directly from the repository: REPLACE WITH THE URL FOR THE MAIN BRANCH

aoi_outline = gpd.read_file('https://raw.githubusercontent.com/nsidc/NSIDC-Data-Cookbook/refs/heads/subsetting-tutorial/example_data/subsetting_tutorial/Himalayas_outline.geojson')
aoi_outline.plot()

Clip the data to the outline with rioxarray

To clip our dataset to our provided outline, we will use rioxarray. rioxarray is an extension of the xarray library that adds geospatial capabilities, allowing Xarray objects to store and work with spatial metadata such as coordinate reference systems (CRS) and spatial dimensions. In other words, it gives Xarray geospatial awareness.

In the following cell we use rioxarray to:
1. Extract the CRS of our dataset, and reproject our region outline to match this projection. Having both datasets in the same CRS ensures everything is spatially aligned before we clip the data.
2. Define the spatial dimensions of our dataset, using the x and y dimensions of ds_subset

# 1. Reproject the outline to the same crs as the data.
aoi_outline = aoi_outline.to_crs(ds_subset.rio.crs)

# 2. Tell rioxarray which dimensions in ds_subset represent the x and y coordinates.
ds_subset = ds_subset.rio.set_spatial_dims(x_dim='x', y_dim='y')

We now have everything in place – rioxarray is aware of the correct spatial dimensions, and our outline is in the same crs as our dataset. We can now put it together to spatially subset the data using rioxarray.clip(). We provide rioxarray.clip() with the desired geometry, the argument all_touched=True to ensure all pixels that overlap with the geometry are included, and the argument drop=True to ensure data outside the clipping area are dropped from the dataset.

ds_clipped = ds_subset.rio.clip(aoi_outline.geometry, all_touched=True, drop=True)
ds_clipped

Our dataset has been successfully clipped to our area of interest! But, you may notice the time coordinate is missing. Same as when we subsetted the dataset by variables, xarray notices the time coordinate isn’t used as a dimension by any of our data variables, so it gets dropped.

We still want the data time holds in our final product, so to fix this we extract the time variable from ds_subset and add it back into ds_clipped using assign_coords().

time_var = ds_subset['time']
ds_clipped = ds_clipped.assign_coords(time=time_var)
ds_clipped.coords

And let’s do a final comparison to see the savings on volume:

ds_clipped_volume = ds_clipped.nbytes
print(f"The beginning volume of the file (uncompressed) was {round((ds_volume / 1e+6), 1)} MBs")
print(f"The volume of the clipped data (uncompressed) is now {round((ds_clipped_volume / 1e+6), 1)} MBs")
print(f"That is a {round(100 - ((ds_clipped_volume / ds_volume)*100),1)} percent reduction in volume for this one file.")

Convert the GeoTransform to GDAL conventions

If we intend to use the clipped dataset in geographic software (such as QGIS), we need to make a small adjustment to its GeoTransform.

A GeoTransform defines the relationship between image pixels and real-world coordinates. As described in the GDAL documentation, it is “an affine transformation from the image coordinate space… to the georeferenced coordinate space.” In simpler terms, it tells GIS software where each pixel belongs on the map.

Let’s first inspect the GeoTransform of the original dataset:

print(ds.rio.transform())

The first six coefficients in the matrix represent:

a – pixel width: 12000.00
b – row rotation: 0.00
c – x-coordinate of the upper-left corner: -3144000.95
d – column rotation: 0.00
e – pixel height: -12000.00
f – y-coordinate of the upper-left corner: -2784000.20

Note that the third (bottom) row of the affine GeoTransform matrix is essentially a placeholder and carries no geographic information itself.

Now compare this with the GeoTransform of the clipped dataset:

print(ds_clipped.rio.transform(recalc=True))

The pixel width and the new origin are what we would expect after clipping. However, notice that the pixel height is now positive (12000) instead of negative (-12000).

This happens because rioxarray derives the GeoTransform directly from the x and y coordinate arrays. In our dataset, both coordinate arrays increase in value, so the resulting transform describes an image whose origin is at the lower-left corner with a positive pixel height.

Many geographic tools, such as GDAL, expect the opposite convention – the origin at the upper-left corner and with a negative pixel height. Meaning, if we exported the dataset at this stage, the data would appear vertically flipped when opened in QGIS.

To make the dataset compatible with GDAL-based software (such as QGIS), we’ll manually construct a new GeoTransform that follows the conventional raster orientation.

Here, we’ll use the Affine constructor to:

  1. Extract the clipped dataset’s GeoTransform, computed by rioxarray.
  2. Calculate the y-coordinate of the upper-left corner.
  3. Create a new GeoTransform with a negative pixel height.
    • With the Affine constructor, the letters a, b, c, etc. correspond to the coefficients labeled above.
    • The pixel width, x-origin, and row/column rotation remains the same, so we pull these as-is from gt.
  4. Write the corrected GeoTransform back to the clipped dataset.
# extract the geotransform of the clipped dataset
gt = ds_clipped.rio.transform(recalc=True)

# calculate y-coordinate of the upper-left corner by multiplying the number of rows by the pixel height and adding this to the y-coordinate of the bottom left corner
y_upper_left = gt.f + ds_clipped.sizes['y'] * gt.e

# create a new geotransform with a negative pixel height and the new y-coordinate of the upper-left corner
new_geotransform = Affine(gt.a, gt.b, gt.c, gt.d, -gt.e, y_upper_left)
print(new_geotransform)

# overwrite the geotransform of the clipped dataset with the new geotransform
ds_clipped = ds_clipped.rio.write_transform(new_geotransform)

Plot the clipped data

Our dataset has been clipped, contains the desired coordinate and data variables, and now holds a GeoTransform that follows GDAL conventions that will georeference the data correctly in QGIS. Let’s plot the same data variable as before to look at our progress.

plot_variable(ds_clipped, 'bc_3d_tot')

5. Export subsetted and clipped dataset as a new NetCDF file

To reiterate a key point: streaming the data directly with earthaccess means no data ever actually touched our local storage or the EC2 disk – everything happened in EC2 memory. At this point in the tutorial, however, our clipped dataset is ready to be written to a new file. We will follow these steps to complete this:

  1. Create a new temp directory that we will store our new file in
  2. Grab the filename of the file we subsetted
  3. Create an output filename by appending the source filename with _subsetted just before the file extension
  4. Write the subsetted dataset to the temp directory
temp_dir = Path("temp")
temp_dir.mkdir(parents=True, exist_ok=True)

filename = results[0]["meta"]["native-id"]
output_filename = Path(filename).with_stem(f"{Path(filename).stem}_subsetted")

ds_clipped.to_netcdf(temp_dir / output_filename, engine='h5netcdf', format='NETCDF4')

6. Move directory containing subsetted files to your local machine

The file HMA2_MATCHA_ATM_20180701T000000Z_V01.0_subsetted.nc has been written to the temp/ folder on our EC2 instance. To move this file to our local machine, we have 2 options: 1. Navigate to the file in the Jupyter Server directory, right click on the file and select ‘Download’ 2. Copy the /temp directory to our local machine using the secure copy (scp) command in the terminal of our local machine.

If the temp/ folder contains many files that we want copied to our local machine (e.g., after looping through the results in section 7 below), then option 2 allows us to do this efficiently. On our local machine, open the terminal and run the command:

scp -i /path/key-pair-name.pem -r ec2-user@instance-public-dns-name:path/to/temp/folder /path/to/local/

Make sure you replace the paths to the pem key, the path to the temp/ folder, the output location, and the name of the EC2 instance with your specific information.

A note on egress costs: When transferring data from an AWS EC2 instance over the public internet (i.e., to your local machine), there may be associated costs with this action known as egress costs. Each month, AWS users can transfer 100GB of data out of AWS EC2 at no cost. As we have seen, we were able to reduce the file to an uncompressed volume of 15.6 MB – well short of the 100GB limit. The cost of transferring data exceeding beyond the free 100GB of data is scaled, starting at $0.09 per GB for the first 10TB. For more information, please see documentation on Amazon EC2 On-Demand Pricing

Lastly, for workflows subsetting large volumes of data that would strain the EC2 instance’s storage, writing the files to an S3 bucket instead can be beneficial. Subsetted data in your S3 bucket can be pulled into the EC2 instance when needed, or downloaded directly to your local machine. Additionally, storing data in S3 is generally cheaper than keeping it on EC2 storage, and will avoid the egress costs noted above until the data is downloaded from S3 to your local machine.

7. Putting it all together

Loop through all results from search_data()

So far, you’ve subsetted and clipped a single file from the search results in section 2. To apply the same steps to every file, loop through results and repeat the subset and clip for each one.
Note - This loop assumes the following objects created earlier in the notebook exist in the current session: - results from earthaccess.search_data() - data_vars, the list of variables to subset on - aoi_outline, the region of interest (and in the same CRS as HMA2_MATCHA)

%%time

total_subset_volume = 0


for result in results:
    
    # grab source filename, and create the output filename
    filename = result["meta"]["native-id"]
    output_filename = Path(filename).with_stem(f"{Path(filename).stem}_subsetted")

    # load the file
    file = earthaccess.open([result])[0]
    
    # open the file using a context manager
    with xr.open_dataset(file, engine='h5netcdf', decode_coords='all', mask_and_scale=False) as ds:

        # subset ds to specified variables
        ds_subset = ds[data_vars]
    
        # subset file to AOI
        ds_subset = ds_subset.rio.set_spatial_dims(x_dim='x', y_dim='y')
        ds_clipped = ds_subset.rio.clip(aoi_outline.geometry, all_touched=True, drop=True)
        # re-insert the time coordinate
        time_var = ds_subset['time']
        ds_clipped = ds_clipped.assign_coords(time=time_var)
        # re-write the geotransform
        gt = ds_clipped.rio.transform(recalc=True)
        y_upper_left = gt.f + ds_clipped.sizes['y'] * gt.e
        new_geotransform = Affine(gt.a, gt.b, gt.c, gt.d, -gt.e, y_upper_left)
        ds_clipped = ds_clipped.rio.write_transform(new_geotransform)

        total_subset_volume += ds_clipped.nbytes
    
        # write the clipped dataset to a new NetCDF file in the temp folder
        ds_clipped.to_netcdf(f"temp/{output_filename}", engine="h5netcdf", format='NETCDF4')
    
    # remove objects from memory and run garbage collection
    del file, ds_subset, ds_clipped
    gc.collect()
    
    print(f"The file {filename} has been processed and a subsetted version saved to the temp folder.")

print("All files have been processed, subsetted, and written to the temp folder.")
print(f"Total uncompressed volume of the {len(results)} subsetted files written to disk: {round((total_subset_volume / 1e+6), 1)} MBs.")

The temp folder now contains 8 subsetted and clipped HMA2_MATCHA NetCDF files. We can now use the scp command from above in our local terminal to move all the files to our local machine for local analysis.