告别Cygwin!用Python+EarthData API搞定MODIS数据自动下载(附完整脚本)
用Python全自动化获取MODIS气溶胶数据的工程实践
当我们需要分析大气污染、气候变化或空气质量时,MODIS气溶胶光学厚度(AOD)数据是不可或缺的基础资料。但传统的数据获取方式往往需要手动操作多个步骤——从浏览器登录、文件筛选到下载管理,整个过程既耗时又容易出错。本文将展示如何用Python构建一个端到端的自动化数据管道,彻底摆脱对图形界面和手工操作的依赖。
1. EarthData账户配置与认证
NASA EarthData是获取MODIS数据的官方门户,但自动化访问需要先完成服务注册和认证配置。不同于传统的手动登录方式,我们将通过程序化手段处理整个认证流程。
首先访问EarthData注册页面创建账户。完成基础注册后,还需要在"Application"页面生成API访问令牌。这个令牌将作为脚本认证的凭证,避免直接存储用户名和密码。
认证配置文件示例(~/.netrc):
machine urs.earthdata.nasa.gov login 您的用户名 password 您的API令牌注意:务必设置文件权限为600,防止凭证信息泄露。在Linux/Mac上执行
chmod 600 ~/.netrc
Python中我们可以使用netrc模块自动处理认证:
import netrc from os.path import expanduser def get_earthdata_creds(): try: secrets = netrc.netrc(file=expanduser("~/.netrc")) return secrets.authenticators("urs.earthdata.nasa.gov") except Exception as e: print(f"认证失败: {str(e)}") return None2. 构建高效的数据查询系统
LAADS DAAC提供了多种查询接口,我们将使用其REST API实现精准的数据筛选。以下关键参数决定了数据获取的精度和效率:
| 参数 | 说明 | 示例值 |
|---|---|---|
| product | 数据产品代码 | MOD04_L2 (气溶胶产品) |
| collection | 数据版本 | 6.1 |
| dateRange | 时间范围 | 2023-01-01,2023-01-31 |
| areaOfInterest | 地理范围 | 110,20,120,30 (经度min,纬度min,经度max,纬度max) |
| dayNightFlag | 日夜标志 | Day |
Python查询函数实现:
import requests from datetime import datetime, timedelta def query_modis_data(start_date, end_date, bbox=None): base_url = "https://ladsweb.modaps.eosdis.nasa.gov/api/v1/content/details" params = { "product": "MOD04_L2", "collection": "6.1", "dateRange": f"{start_date},{end_date}", } if bbox: params["areaOfInterest"] = ",".join(map(str, bbox)) try: response = requests.get(base_url, params=params) response.raise_for_status() return response.json() except requests.exceptions.RequestException as e: print(f"查询失败: {str(e)}") return None这个函数返回的JSON包含了所有匹配文件的基本信息,包括:
- 文件名
- 文件大小
- 下载URL
- 采集时间
- 空间覆盖范围
3. 实现可靠的文件下载管理
直接的大文件下载常会遇到网络中断问题,我们需要实现以下关键功能:
- 断点续传
- 并行下载
- 完整性校验
- 进度显示
增强版下载函数:
import os import hashlib from tqdm import tqdm from concurrent.futures import ThreadPoolExecutor def download_file(url, target_dir, max_retries=3, chunk_size=8192): filename = url.split("/")[-1] filepath = os.path.join(target_dir, filename) tmp_path = f"{filepath}.part" # 检查是否已完整下载 if os.path.exists(filepath): return filepath headers = {} # 处理断点续传 if os.path.exists(tmp_path): downloaded_size = os.path.getsize(tmp_path) headers = {"Range": f"bytes={downloaded_size}-"} for attempt in range(max_retries): try: with requests.get(url, headers=headers, stream=True) as r: r.raise_for_status() total_size = int(r.headers.get("content-length", 0)) + (downloaded_size if headers else 0) mode = "ab" if headers else "wb" with open(tmp_path, mode) as f, tqdm( total=total_size, unit="B", unit_scale=True, desc=filename ) as progress: for chunk in r.iter_content(chunk_size=chunk_size): if chunk: f.write(chunk) progress.update(len(chunk)) # 验证文件完整性 if validate_file(tmp_path): os.rename(tmp_path, filepath) return filepath else: os.remove(tmp_path) continue except Exception as e: print(f"下载失败(尝试 {attempt+1}/{max_retries}): {str(e)}") continue return None def validate_file(filepath): # 实现MD5校验或其他验证逻辑 return True并行下载控制器:
def download_manager(file_list, target_dir, max_workers=4): os.makedirs(target_dir, exist_ok=True) with ThreadPoolExecutor(max_workers=max_workers) as executor: futures = [] for file_info in file_list: url = file_info["downloadsLink"] future = executor.submit(download_file, url, target_dir) futures.append(future) results = [] for future in futures: results.append(future.result()) return [r for r in results if r is not None]4. 数据处理与质量控制
获取原始HDF文件后,我们需要提取有用的科学数据集(SDS)。MODIS气溶胶产品包含多个数据层,最重要的是550nm处的气溶胶光学厚度。
使用PyHDF库读取数据:
from pyhdf.SD import SD, SDC def extract_aod(hdf_path): try: hdf = SD(hdf_path, SDC.READ) sds = hdf.select("AOD_550_Dark_Target_Deep_Blue_Combined") data = sds.get() attrs = sds.attributes() # 应用缩放因子和偏移量 scale = attrs["scale_factor"] offset = attrs["add_offset"] valid_range = attrs["valid_range"] data = data * scale + offset data[data < valid_range[0]] = np.nan data[data > valid_range[1]] = np.nan return { "data": data, "longitude": hdf.select("Longitude").get(), "latitude": hdf.select("Latitude").get(), "qa": hdf.select("AOD_550_Dark_Target_Deep_Blue_Combined_QA").get() } except Exception as e: print(f"读取HDF失败: {str(e)}") return None finally: hdf.end()数据质量控制标志处理:
MODIS数据包含详细的质量控制(QA)信息,我们需要解码这些标志位来筛选可靠数据:
def apply_qa_mask(aod_data, qa_data): # 每个QA字节包含多个标志位 good_quality = np.zeros_like(qa_data, dtype=bool) for i in range(qa_data.shape[0]): for j in range(qa_data.shape[1]): qa = qa_data[i,j] # 检查云污染标志 cloud = (qa >> 0) & 0x03 # 检查气溶胶类型标志 aerosol_type = (qa >> 2) & 0x07 # 检查陆地/水体标志 land_water = (qa >> 5) & 0x03 good_quality[i,j] = (cloud == 0) and (aerosol_type != 7) and (land_water != 3) aod_data[~good_quality] = np.nan return aod_data5. 完整工作流集成
现在我们将各个模块整合成一个完整的自动化管道:
from datetime import datetime, timedelta import json def modis_aod_pipeline(start_date, end_date, bbox=None, output_dir="output"): # 创建输出目录结构 raw_dir = os.path.join(output_dir, "raw_hdf") processed_dir = os.path.join(output_dir, "processed") os.makedirs(raw_dir, exist_ok=True) os.makedirs(processed_dir, exist_ok=True) # 1. 查询可用数据 print("查询MODIS数据目录...") file_list = query_modis_data(start_date, end_date, bbox) if not file_list: print("未找到匹配的数据") return # 2. 下载文件 print(f"开始下载{len(file_list)}个文件...") downloaded_files = download_manager(file_list, raw_dir) # 3. 处理数据 results = [] for hdf_file in downloaded_files: print(f"处理文件: {os.path.basename(hdf_file)}") data = extract_aod(hdf_file) if data: clean_data = apply_qa_mask(data["data"], data["qa"]) results.append({ "aod": clean_data, "lon": data["longitude"], "lat": data["latitude"], "time": extract_time_from_filename(hdf_file) }) # 4. 保存处理结果 output_file = os.path.join(processed_dir, f"aod_{start_date}_{end_date}.nc") save_as_netcdf(results, output_file) print(f"处理完成,结果已保存至{output_file}")定时自动运行设置:
对于长期监测项目,我们可以使用系统定时任务定期执行脚本:
# 每天凌晨2点执行数据获取 0 2 * * * /usr/bin/python3 /path/to/modis_pipeline.py --days 7 --output /data/modis_aod6. 可视化与结果验证
数据处理的最后一步是生成直观的可视化结果,验证数据质量:
import matplotlib.pyplot as plt import cartopy.crs as ccrs def plot_aod(aod_data, lon, lat, title="MODIS AOD 550nm"): plt.figure(figsize=(12, 8)) ax = plt.axes(projection=ccrs.PlateCarree()) ax.coastlines() mesh = ax.pcolormesh(lon, lat, aod_data, vmin=0, vmax=1.5, cmap="Spectral_r", transform=ccrs.PlateCarree()) plt.colorbar(mesh, label="AOD at 550nm", extend="both") ax.set_title(title) plt.tight_layout() return plt在实际项目中,这套自动化系统将数据处理时间从原来的数小时缩短到几分钟,且完全无需人工干预。一个典型的应用场景是空气质量预警系统——当AOD值超过阈值时,系统可以自动触发警报邮件:
def check_aod_threshold(aod_data, threshold=0.8): high_aod = np.nanmean(aod_data) > threshold if high_aod: send_alert_email(aod_data) def send_alert_email(aod_data): import smtplib from email.mime.text import MIMEText avg_aod = np.nanmean(aod_data) msg = MIMEText(f"警告: 检测到高气溶胶浓度,平均AOD值为{avg_aod:.2f}") msg["Subject"] = "空气质量预警" msg["From"] = "aod_monitor@example.com" msg["To"] = "operations@example.com" with smtplib.SMTP("smtp.example.com") as server: server.send_message(msg)通过将这些技术组件有机结合,我们构建了一个健壮、高效的MODIS数据处理管道。这套方案不仅适用于科研场景,也能满足业务化运行的环境监测需求。
