目录
前言
在遥感GIS应用等领域,使用其他作者生产的产品是很常见的,GLASS product是常见的生态、环境、气象领域会应用到的一类遥感产品。这里我们结合python中的arcpy包进行批处理最终获取所需的数据的TIF格式
提示:以下是本篇文章正文内容,下面案例可供参考
一、如何下载GLASS数据
GLASS数据下载的网站:Global LAnd Surface Satellite (GLASS)
或者北师大:http://glass-product.bnu.edu.cn/index.html
这里参考其他文章可以使用Python按年份批量下载方式,这里以下载长江流域的GLASS LAI为例。具体参数需要进行修改。
#导入requests
import os
import requests
#http://www.glass.umd.edu/GPP/MODIS/500m/GLASS_GPP_500M_V60/2000/049/GLASS12E01.V60.A2000049.h25v05.2022059.hdf
base_url = "http://www.glass.umd.edu/LAI/MODIS/500m/"#这个是数据集产品
# 生成年份列表
years = list(range(2001, 2020))#这里我下载的是2000-01-01-2020-12-30的所有数据
# 指定下载数据的保存目录
output_folder = r"H:\GLASS\LAI500m"
#生成数据URL
#该数据集时间分辨率为8天
for year in years:
dirname = r'H:\GLASS\LAI500m'+os.sep+str(year)
if str(year) not in os.listdir(output_folder):
os.mkdir(dirname)
for day in range(1, 362, 8):
urls = []
day_str = str(day).zfill(3)
folder = f"{year}/{str(day).zfill(3)}"
print(folder)
#我的研究区方位为h25v05和h26v05,不同研究区有不同的范围
i = '010'#修改名字
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h25v05.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h26v05.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h27v05.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h28v05.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h25v06.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h26v06.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h27v06.2022{i}.hdf")
urls.append(f"{base_url}{folder}/GLASS01E01.V60.A{year}{day_str}.h28v06.2022{i}.hdf")
# 下载数据到指定的文件夹
for url in urls:
response = requests.get(url)
if response.status_code == 200:
filename = url.split("/")[-1]
filepath = output_folder+os.sep+str(year)+os.sep+filename
with open(filepath, "wb") as f:
f.write(response.content)
print(f"Downloaded {filename}")
else:
print(f"Failed to download {url}")
另外网络上还有使用IDM下载器进行处理的操作具体可以参考:小工具 | 网站数据抓取(以GLASS数据为例)-腾讯云开发者社区-腾讯云
二、下载后的GLASS数据如何处理——以LAI为例
我们的最终目的是得到研究区的逐月GLASS数据(TIF格式)
1.影像格式转换
代码如下(示例):将文件从hdf格式转为TIF格式
# -*- coding: UTF-8 -*-
import arcpy
import os
for time in range(2018,2022):
#批量读取一级目录下的子文件
rootdir = r"H:\GLASS\LAI500m" + os.sep +str(time)
# rootfiles = os.listdir(rootdir)
hdf_file_name_list = os.listdir(rootdir)
for hdf_file in hdf_file_name_list:
if hdf_file.endswith(".hdf"):
tif_file_name = hdf_file[0:30] + ".tif"
print (tif_file_name)
arcpy.ExtractSubDataset_management(rootdir + os.sep + hdf_file, rootdir + os.sep + tif_file_name)
print (rootdir + os.sep + tif_file_name + "完成")
2.影像镶嵌
将下载的每一期的影像拼接起来得到整幅涵盖研究区范围的无缝影像,当然后续还需要进行影像裁剪操作以匹配研究区范围。
# -*- coding: UTF-8 -*-
import arcpy
import glob
import os
import datetime
arcpy.CheckOutExtension("Spatial") # 检查是否注册该模块
def batchmosaic(time):
# 输入路径,这里选中上一步解压后保存解压数据的文件夹
inws = r"H:\GLASS\LAI500m"+os.sep + str(time)
# 输出路径
outws = r"H:\GLASS\LAI_mosaic" + os.sep + str(time)
for day in range(1, 362, 8):
urls = []
day_str = str(day).zfill(3)
# 不同研究区有不同的范围需要自行设置,数量多的话可进行for循环
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h25v05.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h26v05.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h27v05.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h28v05.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h25v06.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h26v06.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h27v06.tif'%(str(time),day_str))
urls.append(inws+os.sep+'GLASS01E01.V60.A%s%s.h28v06.tif'%(str(time),day_str))
# 以下逻辑为:
# 获得子目录全路径,并进入子目录
base = urls[0]
print(base)
out_coor_system = arcpy.Describe(base).spatialReference # 坐标系统
dataType = arcpy.Describe(base).DataType # 数据类型
piexl_type = arcpy.Describe(base).pixelType # 像素类型
cellwidth = arcpy.Describe(base).meanCellWidth # 像元宽度
bandcount = arcpy.Describe(base).bandCount # 获取bandCount
mosaic_ras = [] # 待拼接影像的文件名,且中间以;隔开,例如:a.tif;b.tif;c.tif
for ras in urls:
mosaic_ras.append(ras)
mosaic_ras = ";".join(urls) # 字符串拼接
outFolder = outws+os.sep+day_str
isExists = os.path.exists(outFolder)
if not isExists:
os.makedirs(outFolder)
nameT = '%s_%s.tif'%(str(time),day_str)
# 执行拼接操作
arcpy.MosaicToNewRaster_management(mosaic_ras, outFolder, str(nameT), out_coor_system,
"32_BIT_FLOAT",
cellwidth, bandcount, "LAST", "FIRST")
print(os.path.basename(nameT) + " ---- OK! ---- ")
print(inws + "has Finish!")
print("-----------------------------------------------------------")
for time in range(2018,2022):
batchmosaic(time)
3.影像裁剪与投影
代码如下(示例):记得修改掩膜、输入、输出路径以及坐标系,这里使用的Web Mercator投影坐标系,
import arcpy
import os
Coordinate_System = "PROJCS['WGS_1984_Web_Mercator_Auxiliary_Sphere',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Mercator_Auxiliary_Sphere'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',0.0],PARAMETER['Standard_Parallel_1',0.0],PARAMETER['Auxiliary_Sphere_Type',0.0],UNIT['Meter',1.0]]"
inputdir=r"H:\GLASS\LAI_mosaic"
output_path = r"H:\GLASS\LAI500m\mosaic_clip"
mask = "H:\database\Yangtze.shp"
arcpy.env.parallelProcessingFactor = '0'
for time in range(2018,2022):
seconddir = inputdir + os.sep + str(time)
thirddirs = os.listdir(seconddir)
for thirddir in thirddirs:
wholedir = seconddir + os.sep + thirddir
print wholedir
arcpy.env.workspace = wholedir
rasterList = arcpy.ListRasters("*", "TIF")
print rasterList
for raster in rasterList:
outname = output_path + os.sep + str(time)+raster[5:8] + "Web84" + ".tif"
if os.path.exists(outname):
print (outname + "have exist")
else:
arcpy.ProjectRaster_management(raster, outname, Coordinate_System, "NEAREST", "#", "#", "#", "#")
print (outname + "have been done")
out = output_path + os.sep + str(time)+raster[5:8] + ".tif"
arcpy.Clip_management(outname, "#", out, mask, "-999", "ClippingGeometry","NO_MAINTAIN_EXTENT")
print out + " clip has done!"
4.影像最大值合成
由于得到的GLASS影像是以8天为时间间隔,所以很多模型进行分析是以月为尺度单位,所以需要对每一月的数据进行合成,这一步就得到的了研究区月度的GLASS产品数据的TIF格式文件了。
具体的操作可以参考之前的文章:
北师大GLASS LAI数据的处理流程及使用北师大GLASS LAI数据进行月度数据合成(最大值)_JianShu Wang的博客-CSDN博客
5.影像比例系数缩放
由于下载得到的数据一般都有一个比例系数的缩放,也就是产品下载后产品的值乘上比例系数才是真实的数值,可以进行栅格计算得到真实结果,之后产品便可以进行分析和使用了。
总结
这里以GLASS LAI500m2018-2021年为例,在马里兰大学的官网进行了下载与数据处理,最终得到可进行使用的TIF格式产品数据的整个预处理路线。如有问题欢迎提出!
注意-如果代码跑不通极大可能是文件路径或文件名设置的问题,整体的逻辑是没有问题的。