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

目录

前言

一、如何下载GLASS数据

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

1.影像格式转换 

2.影像镶嵌

3.影像裁剪与投影

4.影像最大值合成

5.影像比例系数缩放

总结


前言

在遥感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格式产品数据的整个预处理路线。如有问题欢迎提出!

注意-如果代码跑不通极大可能是文件路径或文件名设置的问题,整体的逻辑是没有问题的。

  • 8
    点赞
  • 58
    收藏
    觉得还不错? 一键收藏
  • 17
    评论
Glass数据集是一个用于分类的公开数据集,共有214个样本,每个样本有10个特征。该数据集用于根据样本的特征预测玻璃的类型,有7种可能的玻璃类型。 要对Glass数据集进行分类,可以使用各种机器学习算法和技术,如决策树、支持向量机、随机森林等。以下是一种基本的分类方法: 首先,需要将数据集划分为训练集和测试集,以便评估模型的性能。可以将数据集的70-80%用作训练集,剩余的30-20%用作测试集。 然后,选择适当的机器学习算法。一种常见的分类算法是决策树。决策树通过使用特征值逐步划分数据,构建出一颗树状结构,从而对新样本进行分类。 接下来,使用训练集训练决策树模型。训练过程中,根据样本的特征值进行划分,并根据标签值对每个划分进行评估,选择最佳划分。 训练完成后,使用测试集评估模型的性能。将测试集的样本输入到训练好的决策树模型中,观察模型对样本的分类结果,并与真实标签进行比较,计算准确率、召回率等指标来评估模型的性能。 如果模型的性能不尽如人意,可以尝试使用其他算法或者进行参数调整来提高分类的准确性。可以尝试使用支持向量机、随机森林等算法,并进行交叉验证来选择最优的参数。 最后,可以使用整个数据集进行模型训练,以获取最终的分类模型。使用该模型可以对新的未知样本进行分类预测,从而实现对Glass数据集的分类任务。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 17
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值