How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?

Published: October 19, 2024

Updated: October 19, 2024

Tags: MODIS; NASA; Python;

DMCA.com Protection Status

Introduction

This article demonstrates how to download and plot NASA MODIS L3 (Level-3) products over polar regions. We'll cover downloading the data from NASA's LAADS DAAC, extracting key scientific datasets, and creating visualizations using Python tools like Cartopy, Matplotlib, and PyHDF.

Import Python Libraries

First, we need to import the necessary Python libraries to download, read, and visualize the data.

1
2
3
4
5
6
7
8
import urllib.request, json
import os
import netCDF4
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from datetime import date
from pyhdf.SD import SD, SDC

Cartopy is particularly useful for mapping geographic data, and pyhdf enables reading of MODIS HDF files.

Download Data from NASA LAADS DAAC

To get NASA MODIS data, we use the LAADS DAAC API. Here, we download monthly MODIS Level-3 cloud optical thickness data for August 2019 from the Aqua satellite. Ensure you have an Earthdata account and a valid authentication token.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
output_dir = 'Datasets'
token = 'your_token_here'
year, month = 2019, 8
platform, instrument, product_category, product = 'AQUA', 'MODIS', 61, 'MYD08_M3'

# Calculate day count for accessing data by date
d0, d1 = date(year, 1, 1), date(year, month, 1)
count_of_day = (d1 - d0).days + 1

# Construct the LAADS DAAC URL for the product
ladsweb_url = f'https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/{product_category}/{product}/{year}/{count_of_day}.json'

Next, download and save the MODIS HDF file.

1
2
3
4
5
6
7
8
with urllib.request.urlopen(ladsweb_url) as url:
    data = json.loads(url.read().decode())

for file in data['content']:
    file_name = file['name']
    target_dir = f'{output_dir}/{platform}/{instrument}/{product}/{year}/{year}_{month}_01/{file_name}'
    if not os.path.exists(target_dir):
        urllib.request.urlretrieve(f'https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/{product_category}/{product}/{year}/{count_of_day}/{file_name}', target_dir)

Reading MODIS Data Using PyHDF

Once the data is downloaded, we can use the pyhdf library to explore and extract datasets from the HDF file.

1
2
3
file_name = 'Datasets/AQUA/MODIS/MYD08_M3/2019/2019_08_01/MYD08_M3.A2019213.061.2019244183905.hdf'
file = SD(file_name, SDC.READ)
file.info()

This command outputs the number of datasets (1144) in the file. Let’s search for specific datasets related to cloud optical thickness.

1
2
3
datasets_dic = file.datasets()
words = ['Cloud', 'Optical', 'Thickness', 'Liquid']
matches = [dataset for dataset in datasets_dic.keys() if all(word in dataset for word in words)]

The matches contain cloud optical thickness datasets like Cloud_Optical_Thickness_Liquid_Mean_Mean, which we will use for plotting.

Extracting Cloud Optical Thickness Data

To extract the data and apply the scale factor:

1
2
3
4
5
6
7
sds_obj = file.select('Cloud_Optical_Thickness_Liquid_Mean_Mean')
data = sds_obj.get()
attributes = sds_obj.attributes()
scale_factor = attributes['scale_factor']

# Apply scale factor to correct the data
data = data * scale_factor

The resulting data array contains cloud optical thickness values that range from 0 to 80 for valid measurements, while invalid measurements are represented by the fill value (-9999).

Plotting MODIS Data

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
fig = plt.figure(num=None, figsize=(8, 6), dpi=80, edgecolor='k')

ax = plt.axes(projection=ccrs.PlateCarree())

im = ax.imshow(data, extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='gist_earth',vmin=0,vmax=80)

ax.coastlines()

plt.title("Cloud_Optical_Thickness_Liquid_Mean_Mean", fontsize=12)

plt.colorbar(im,fraction=0.02)

#plt.savefig("cartopy_heatmap_02.png", bbox_inches='tight', dpi=200)

How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?
How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?

Plotting MODIS Data Over the Arctic

Example 1: Basic Arctic Plot

Here, we visualize the data over the entire Arctic region using Cartopy's NorthPolarStereo projection.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
fig = plt.figure(num=None, figsize=(8, 6), dpi=80, edgecolor='k')

ax = plt.axes(projection=ccrs.NorthPolarStereo())

ax.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())

im = ax.imshow(data, extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='gist_earth', alpha=0.6,vmin=0,vmax=20)

ax.coastlines()

plt.colorbar(im)

plt.title('Cloud_Optical_Thickness_Liquid_Mean_Mean', 
          fontsize=10)

# Save the figure as a high-resolution PNG file and display it
plt.savefig("cartopy_arctica_03.png", bbox_inches='tight', dpi=200)

plt.show()

How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?
How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?

Example 2: Enhanced Polar Plot with Circular Cropping

For a more refined visualization, we can crop the map into a circular region and display only the Arctic region.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
# Create figure and axis using North Polar Stereographic projection
fig = plt.figure(figsize=[10, 5])
ax1 = plt.subplot(1, 2, 1, projection=ccrs.NorthPolarStereo())

# Set map extent: focusing on latitudes between 60 and 90 degrees (the Arctic region)
ax1.set_extent([-180, 180, 60, 90], ccrs.PlateCarree())

# Add coastlines for geographic context
ax1.coastlines()

# Add gridlines to show latitude and longitude markers
ax1.gridlines()

# Define a circular boundary to crop the map within a circular region
theta = np.linspace(0, 2 * np.pi, 100)
center, radius = [0.5, 0.5], 0.5
verts = np.vstack([np.sin(theta), np.cos(theta)]).T
circle = mpath.Path(verts * radius + center)

# Apply the circular boundary to the polar map to focus on the Arctic region
ax1.set_boundary(circle, transform=ax1.transAxes)

# Add title
plt.title('Cloud_Optical_Thickness_Liquid_Mean_Mean', 
          fontsize=10)

im = ax1.imshow(data, extent=[-180,180,-90,90], transform=ccrs.PlateCarree(),cmap='gist_earth', alpha=0.6,vmin=0,vmax=20)

# Save the figure as a high-resolution PNG file and display it
plt.savefig("cartopy_arctica_04.png", bbox_inches='tight', dpi=200)
plt.show()

How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?
How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?

Plotting MODIS Data Over the Antarctica

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
# Create the figure
fig = plt.figure(num=None, figsize=(8, 6), dpi=80, edgecolor='k')

# Set up the South Polar Stereographic projection
ax = plt.axes(projection=ccrs.SouthPolarStereo())

# Set the extent to focus on Antarctica
ax.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())

# Plot the data with imshow using the 'gist_earth' colormap and a specific range for values (vmin, vmax)
im = ax.imshow(data, extent=[-180, 180, -90, 90], transform=ccrs.PlateCarree(),
               cmap='gist_earth', alpha=0.6, vmin=0, vmax=80)

# Add coastlines for better geographical context
ax.coastlines()

# Add color bar for the plot
plt.colorbar(im)

# Add the title for the plot
plt.title("Cloud_Optical_Thickness_Liquid_Mean_Mean", fontsize=12)

# Display the plot
plt.savefig("Cloud_Optical_Thickness_Liquid_Mean_Mean_2019_08_fig_04.png", bbox_inches='tight', dpi=200)
plt.show()

How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?
How to Plot NASA MODIS L3 Products Over Polar Regions using Python ?

References

Links Site
colormaps matplotlib.org
L3 MODIS Monthly Imagery atmosphere-imager.gsfc.nasa.gov
Image

of