MODIS SMI数据批处理

概述

最近在做海洋遥感课程的大作业：分析西北太平洋台风对Chla、SST、SSH的影响。关于如何分析台风的影响是不知道了，只打算把台风周期内每天的Chla和SST出图，然后定性描述下Chla和SST是如何变化的，暂且如此。
文章主要记录从数据下载到出图的一个完整流程。

数据批量下载：

数据使用NASA Level-3 Standard Mapped Image Products，数据说明详见用户手册
将下载链接保存到本地：

# Python 3.6.6
#保存下载链接
import os,requests,re
from bs4 import BeautifulSoup as bs

URL = 'https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/Mapped/Daily/9km/chlor_a/2016/'  #链接地址
DIR = 'E:\data\chla'  #保存路径
FILENAME = 'chla_2016.txt' #保存文件
def WriteURL(url,filename):
	resp = requests.get(url)
	pattern = r'<td><a href=\'(.*?)\'>'
	down_url_list = re.findall(pattern,resp.text)
	with open(filename,'w') as file:
		for down_url in down_url_list:
			text = down_url + '\n'
			print(text)
			file.write(text)
			file.flush()
if __name__=='__main__':
	os.chdir(DIR) #
	WriteURL(URL,FILENAME)

该数据下载的时候无需登录，可以直接下载：

import wget,os,glob
DIR = 'E:\data\chla' #路径

def getUrls():
	urls = []
	UrlFile = glob.glob('./*.txt')[0]
	with open(UrlFile,'r') as f:
		for line in f:
			url = line.replace('\n','')
			urls.append(url)
	return urls

if __name__=='__main__':
	os.chdir(DIR)
	for url in getUrls():
		name = url.split('/')[-1]
		if os.path.exists(name):
			continue
		wget.download(url)

同理可将2016年SST的数据下载下来。

HDF数据转为GeoTiff

下载到的数据后缀名为‘.nc，按理说应该是NetCDF格式的文件，然而在用户手册中提到文件格式为HDF5，但在ArcGIS中可以NetCDF格式直接打开，又可以HDF5格式在Python中操作，很是有趣。这里以SST为例，在Python中将数据转为GeoTIff格式：

# Python 3.6.6
from osgeo import gdal
from osgeo import osr
import numpy as np
import os,sys,glob,h5py
DIR = 'e:\ocean\MODIS_SST'

def NcToRaster(file):
	# read in data
	with h5py.File(file,'r') as f:
		sst =  f['sst'][:]*0.005 #scale_factor=0.005,add_offset=0
		#真实数据=原始数据*scale_factor+add_offset
		#不同产品scale&offset可能不一致，可以在属性中查看。
		lon = f['lon'][:]
		lat = f['lat'][:]	
	# set geotransform
	nx = sst.shape[0]
	ny = sst.shape[1]
	xmin, ymin, xmax, ymax = [lon.min(), lat.min(), lon.max(), lat.max()] 
	xres = (xmax - xmin)/float(ny) 
	yres = (ymax - ymin)/float(nx) 
	geotransform = (xmin, xres, 0, ymax, 0, -yres) 	
	# create raster file
	filename = file+'.tif'
	dst_ds = gdal.GetDriverByName('GTiff').Create(filename, ny, nx, 1, gdal.GDT_Float32)这里需要注意数据类型应该和原始数据一致
	dst_ds.SetGeoTransform(geotransform) # specify coords 
	srs = osr.SpatialReference() # establish encoding 
	srs.ImportFromEPSG(4326)   # WGS84 lat/long  https://epsg.io/4326
	dst_ds.SetProjection(srs.ExportToWkt()) # export coords to file 
	dst_ds.GetRasterBand(1).WriteArray(sst) #write sst to raster
	dst_ds.FlushCache() #write to disk
	dst_ds = None #save,close
if __name__=='__main__':
	os.chdir(DIR)
	for file in glob.glob('.\*.nc'):
		NcToRaster(file)

＃批量裁剪

# Python 2.7.12
import arcpy
import glob
import os

inws = r'E:\ocean\MODIS_SST' #输入路径
outws = r'E:\ocean\SST_clip' #输出路径
mask = r'E:\ocean\mask\ResearchArea.shp' #裁剪区域
raster = glob.glob(os.path.join(inws,'*.tif')) #需要批量裁剪的文件
for ras in raster:
    outname = os.path.join(outws,os.path.basename(ras).split('.')[0])+'_SST_clip.tif'  #保存文件名
    arcpy.Clip_management(ras,'#',outname,mask,'-163.835','NONE','NO_MAINTAIN_EXTENT') #执行裁剪

批量出图

# Python 2.7.12
import arcpy,os,glob

inws = r'E:\ocean\chla_clip' #输入路径
outws = r'E:\ocean\mapping\chla' #输出路径

tifs = glob.glob(os.path.join(inws,'*.tif')) #需要出图的文件
print(tifs)
mxd = arcpy.mapping.MapDocument(r'E:\ocean\mapping.mxd') #先在mxd设置同意的模板
df= arcpy.mapping.ListDataFrames(mxd)[0] #获取第一个DataFrame
lyr = arcpy.mapping.ListLayers(mxd)[4]   #保留前四个图层不变，更换第五个图层
for tif in tifs:
    paths = os.path.split(tif) #获取文件路径及文件名
    outname = os.path.join(outws,paths[1].replace('_clip.tif','.png')) 输出图名
    lyr.replaceDataSource(paths[0],'NONE',paths[1]) #更换数据源
    arcpy.mapping.ExportToPNG(mxd,outname) #出图
print("OK")

出图

其中台风路径数据下载自中国气象局热带气旋资料中心。