记录GLASS数据从下载到处理——获取研究区逐月GLASS产品

JianShu Wang

已于 2023-06-27 14:32:26 修改

阅读量5k

点赞数 9

文章标签： python arcgis 经验分享学习

于 2023-06-27 14:31:46 首次发布

本文链接：https://blog.csdn.net/DaShuJuZhang/article/details/131415006

版权

前言

一、如何下载GLASS数据

二、下载后的GLASS数据如何处理——以LAI为例

前言

在遥感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!"