实验下载了大量的modis的NDVI数据,由于数据是HDF格式,无法直接在arcgis里镶嵌拼接,同时因为是月数据,工作量很大。基于arcgis自带的arcpy
也可以直接使用arcgis对数据格式进行转换,使用镶嵌工具进行镶嵌:
再使用镶嵌至新栅格进行镶嵌,
也可以使用代码
代码也是参考的别人的,能保证格式转换的方法没问题,下面的影像拼接方法,我没有去验证,后来是使用的其他方法进行的镶嵌。
从这个博文中参考来的,中间还要生成天数列表等,自己没有弄
http://t.csdn.cn/sQ6Uo
首先对数据集进行格式转换:
# coding=utf-8
import os
import arcpy
from arcpy import env
sourceDir = 'F:/NPP2020' # 输入
targetDir = 'F:/NPP2020TIF' # 输出
#检查Spatial —ArcGIS Spatial Analyst 扩展模块是否许可
arcpy.CheckOutExtension("Spatial")
#设置工作空间
env.workspace = sourceDir
arcpy.env.scratchWorkspace = sourceDir
#读取栅格列表
hdfList = arcpy.ListRasters('*', 'hdf')
# for 循环
for hdf in hdfList:
# print hdf #打印文件名用于检查
# 替换文件名后缀
NDVIName = os.path.basename(hdf).replace('hdf', 'tif')
outname = targetDir + '\\' + NDVIName
print(outname) #打印文件名用于检查,python2或3的打印方式不同
arcpy.ExtractSubDataset_management(hdf, outname)
# 设置为0,获取NDVI,设置为1获取EVI层,等等
之后根据每个栅格数据的时间对影像进行拼接:
# coding=utf-8
# -*- coding:gbk-*-
import arcpy
import glob
import os
from arcpy import env
import pandas as pd
arcpy.CheckOutExtension("Spatial")
# 设置工作空间,即上一步tif影像的保存目录
env.workspace = r'D:\MOD13A1TIF' # 输入路径
#输入路径
inws = r'D:\MOD13A1TIF' # 输入路径
##inws=inws.encode('utf-8')
outws = r'D:\MOD13A1UNION' # 输出路径
# 投影文件,可以在arcgis里导出
coordinate = r"D:\WGS84UTMZONE50N.prj"
#days_inws = r'E:\tianshuliebiao'
# 判断日期函数,由于MODIS数据文件是以年内天数命名,为了将其转为日期,并给镶嵌后的文件重新命名,如2001001变为2001_01_01
def out_date(year, day):
fir_day = datetime.datetime(year, 1, 1)
zone = datetime.timedelta(days=day - 1)
return datetime.datetime.strftime(fir_day + zone, "%Y_%m_%d")
# 利用glob包,将inws下的所有tif文件读存放到rasters中
rasters_file = glob.glob(os.path.join(inws, "*.tif"))
##csv_file = glob.glob(os.path.join(days_inws, "*.csv"))
for year in raster_file:
year_csv = os.path.basename(year).split('.')[0]
data = pd.read_csv(year, encoding='gbk')
days = list(data['days'])
for d in days:
raster_list = []
day = int(str(d))
print day, year_csv
for ras in arcpy.ListFiles("*.tif"):
# 循环rasters中的所有影像,进行按掩模提取操作
info = os.path.basename(ras).split('.')[1]
info = info.replace('A', '')
day_info = int(info[4:7])
if str(info[:4]) == str(year_csv) and day == day_info:
raster_list.append(ras)
year__info__shift = info[:4]
days__info__shift = info[4:7]
# 生成文件名日期
finame__1 = out_date(int(year__info__shift), int(days__info__shift))
print( finame__1)
## #完整文件名
outname = finame__1 + '_MOD13A3.tif'
print(outname)
try:
#采用最大值合成法,镶嵌至新栅格
arcpy.MosaicToNewRaster_management(raster_list, outws, outname, coordinate, "32_BIT_UNSIGNED",None, "1",
"MAXIMUM", "FIRST")
except:
print "Mosaic To New Raster example failed."
print arcpy.GetMessages()
print 'mosaic finish'
也可以使用其他工具,比如extract subdataset dengdeng gongju