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) |
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() |
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() |
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() |
References
Links | Site |
---|---|
colormaps | matplotlib.org |
L3 MODIS Monthly Imagery | atmosphere-imager.gsfc.nasa.gov |