前言

最近需要下载气象数据——CHIRPS，借助之前学的批量下载哨兵轨道数据和批量解压，特此写一个脚本提取所需要的数据

一、CHIRPS是什么？

CHIRPS（Climate Hazards Group InfraRed Precipitation with Station data）是一个全球范围的降雨数据集，结合了0.05°分辨率的卫星图像和原位站数据，形成网格化的降雨时间序列，用于趋势分析和季节性干旱监测。数据从1981年开始，可以分为每日/每月/每年，其中又可以细分为每若干月、每若干年等等，详细介绍网站：https://www.chc.ucsb.edu/data/chirps
数据下载网站：https://data.chc.ucsb.edu/products/CHIRPS-2.0/

数据为全球范围数据，坐标系为WGS:84

二、实现步骤

1.下载数据

CHIRPS数据的下载与哨兵轨道数据下载方式略有差别，为了加快下载速度我调用了IDM进行数据的下载，这里需要自行安装一个IDM应用程序，正常来说IDM是有使用期限的，分享一个绿色版的链接，如有需要请自己网盘下载。

idm官方命令参数介绍，如下：

idman /s
或 idman /d URL [/p 本地_路径] [/f 本地_文件_名] [/q] [/h] [/n] [/a]参数：
/d URL - 下载一个文件，等等。
IDMan.exe /d "http://www.internetdownloadmanager.com/path/File Name.zip"
/s - 开始任务调度里的队列
/p 本地_路径 - 定义要保存的文件放在哪个本地路径
/f 本地local_文件_名 - 定义要保存的文件到本地的文件名
/q - IDM 将在成功下载之后退出。这个参数只为第一个副本工作
/h - IDM 将在成功下载之后挂起您的连接
/n - 当不要 IDM 询问任何问题时启用安静模式
/a - 添加一个指定的文件 用 /d 到下载队列，但是不要开始下载参数 /a, /h, /n, /q, /f 本地_文件_名， /p 本地_路径 工作只在您指定文件下载 /d URL

实现代码如下：

    IDM = "E:\IDM\IDM\IDMan.exe"command = ' '.join([IDM, '/d', DownUrl, '/p', DownPath, '/f', FileName, '/a'])os.system(command)call([IDM, '/s'])

需要改的地方为IDM这个变量，这个变量后面的路径为你安装该程序的位置，command是一种类似命令行的调用，call代表开始下载，代码无误应该是下图所示的样子，然后打开IDM可以看到下载队列了。

2.解压缩

在写这部分代码时自己发现一个问题，那就是python调用IDM进行下载后无法实时返回下载情况，1中的代码主要实现：添加下载链接、开始下载，没有返回的信息，而对文件进行解压时就必然要对文件路径进行遍历。因此在IDM执行下载与文件解压之间需要有一段缓冲时间，这里我用通过需要下载文件的个数与实际已经下载好的文件个数进行对比，返回布尔值来决定是否进行下一步的操作，代码如下（示例）：

def Judge_check(array,path):file_list = os.listdir(path)file_array = []for file in file_list:if os.path.splitext(file)[-1] == '.gz':file_array.append(file) #统计已 .gz结尾的文件'=if len(array) == len(file_array):return  Trueelse:nowtime = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime())print("文件尚未下载完全,缺失%d个,120S后重新读取，%s"%((len(array)-len(file_array)),nowtime))return Falsedef Check_file(array,path):file_list = os.listdir(path)file_array = []for file in file_list:if os.path.splitext(file)[-1] == '.gz':file_array.append(file) #统计 .gz结尾的文件if file in array:print("%s完成下载"%file)if len(file_array) == len(array):print("所需时间段内气象数据下载完毕,共计%d个文件"%len(file_list))#exit(
......for i in range(1000):Judge = Judge_check(Chirps_Time_Info,out_path)if Judge == True :Check_file(Chirps_Time_Info,out_path)breakelse:time.sleep(120)#120S后重新读取文件夹下文件信息continue

这里用到一个for循环，通过Judge函数判断文件路径下有多少个已经下载的文件夹，如果没满足要求就返回False并等待120秒后跳出此处循环（continue）再次进行Judge函数判断，直至二者数量对的上返回True，并执行Check函数，Check函数主要统计以.gz结尾的文件信息

然后再执行解压缩操作，具体通过调用gzip库以及shutil库中的函数进行，这里我卡住挺久，具体原因就是IDM下载时会新建一个文件夹，文件夹下才是我们需要的数据，因此路径不是os.path.join(path,file)而是os.path.join(path,file,file)，如果用os.path.join(path,file)那么在open(gz_out_path, “wb+”)就会报错，具体原因就是只有文件才能写入，文件夹是不能写入的

for file in file_list :#!!只有文件才能写入if os.path.splitext(file)[-1] == '.gz':gz_path = os.path.join(path,file,file)f_file = file.replace(".gz", "")# 获取文件的名称，去掉后缀名gz_out_path = os.path.join(path,'Chirps',f_file)#读取解压后的文件，并写入去掉后缀名的同名文件if os.path.exists(gz_out_path) == False:with gzip.open(gz_path,"rb") as f_in:with open(gz_out_path, "wb+") as f_out :shutil.copyfileobj(f_in,f_out)print("%s解压完成"%file)

代码无误后应该显示为下图所示：

3.批量裁剪

这步主要依靠Gdal库和shape库实现，还有一些小的问题，但是不影响最终结果，后续改进的话可以采用更为方便的Geopandas库实现,原理应该是类似arcgis的掩膜提取功。
教训：如果不加cropToCutline=True，那么是矢量裁剪后的tif影像会有偏移…困扰了自己一晚上

ds = gdal.Warp(out_raster,# 裁剪图像保存完整路径（包括文件名）in_raster,# 待裁剪的影像format = 'GTiff',# 保存图像的格式cutlineDSName = in_shape, # 矢量文件的完整路径cutlineWhere=in_shape,# optionally you can filter your cutline (shapefile) based on attribute valuescropToCutline=True, # 将目标图像的范围指定为cutline 矢量图像的范围dstSRS='EPSG:4326', # 参考：WGS84#cutlineWhere="FIELD = 'whatever'", # clip specific feature#以shp的范围输出outputBounds=r.bbox,copyMetadata=True,dstNodata = 0)              # select the no data value you likeds=None     #do other stuff with ds object, it is your cropped dataset. in this case we only close t

小问题就是他会报错，显示:SQL Expression Parsing Error: syntax error, unexpected $undefined, expecting end of string,但是不影响数据，注释掉cutlineWhere=in_shape即可不报错（我也忘记为什么之前加上了）

三、完整代码如下

import urllib.request
from bs4 import BeautifulSoup
import re
import os
import requests
from subprocess import call
import gzip
import shutil #复制文件
import time
import datetime
import shapefile
from osgeo import gdal
import sys""" Notice
# URL https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/chirps-v2.0.1981.02.tif.gz
数据是WGS84
"""def download(dest_dir, url):print(url)print(dest_dir)headers={"Accept": "text/html,application/xhtml+xml,application/xml;q=0.9,image/webp,image/apng,*/*;q=0.8,application/signed-exchange;v=b3;q=0.9","Accept-Encoding": "gzip, deflate, br","Accept-Language": "zh-CN,zh;q=0.9,en;q=0.8,en-GB;q=0.7,en-US;q=0.6","Connection": "max-age=0","Host": "data.chc.ucsb.edu","sec-ch-ua": "\" Not A;Brand\";v=\"99\", \"Chromium\";v=\"100\", \"Google Chrome\";v=\"100\"","sec-ch-ua-mobile": "?0","sec-ch-ua-platform": "\"Windows \"","Sec-Fetch-Dest": "document","Sec-Fetch-Mode": "navigate","Sec-Fetch-Site": "none","Sec-Fetch-User": "?1","Upgrade-Insecure-Requests": "1","User-Agent": "Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/96.0.4664.110 Safari/537.36",}try:request = urllib.request.Request(url,headers=headers)response = urllib.request.urlopen(request)urllib.request.urlretrieve(url,dest_dir)except:error_url.append(url)print("\tError retrieving the URL:", dest_dir)else: # 没有异常if url in error_url: #在错误列表里error_url.remove(url)def IDMdownload(DownUrl, DownPath, FileName):IDM = "E:\IDM\IDM\IDMan.exe"command = ' '.join([IDM, '/d', DownUrl, '/p', DownPath, '/f', FileName, '/a'])os.system(command)call([IDM, '/s'])def Judge_check(array,path):file_list = os.listdir(path)file_array = []for file in file_list:if os.path.splitext(file)[-1] == '.gz':file_array.append(file) #统计已 .gz结尾的文件'=if len(array) == len(file_array):return  Trueelse:nowtime = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime())print("文件尚未下载完全,缺失%d个,120S后重新读取，%s"%((len(array)-len(file_array)),nowtime))return Falsedef Check_file(array,path):file_list = os.listdir(path)file_array = []for file in file_list:if os.path.splitext(file)[-1] == '.gz':file_array.append(file) #统计 .gz结尾的文件if file in array:print("%s完成下载"%file)if len(file_array) == len(array):print("所需时间段内气象数据下载完毕,共计%d个文件"%len(file_list))#exit(def ungzip(path):file_list = os.listdir(path)store_path = os.path.join(path + './Chirps')if os.path.exists(store_path) == False: #布尔逻辑词不能用字符os.mkdir(path + './Chirps')for file in file_list :#!!只有文件才能写入，文件夹是不能写入的if os.path.splitext(file)[-1] == '.gz':gz_path = os.path.join(path,file,file)f_file = file.replace(".gz", "")# 获取文件的名称，去掉后缀名gz_out_path = os.path.join(path,'Chirps',f_file)#读取解压后的文件，并写入去掉后缀名的同名文件if os.path.exists(gz_out_path) == False:with gzip.open(gz_path,"rb") as f_in:with open(gz_out_path, "wb+") as f_out :shutil.copyfileobj(f_in,f_out)print("%s解压完成"%file)else:print("%s已解压完成，无需再次解压"%file)continuef_out.close()f_in.close()else:continue #跳过循环def warpclip(in_shape,in_dir,out_dir):in_shape = in_shape #shp路径in_dir = in_dir # 输入栅格图像地址out_dir = out_dir # 输出裁剪影像地址file_list = os.listdir(in_dir) #获取输入图像文件夹的地址r = shapefile.Reader(in_shape,encoding='gb18030') #获取shp的范围并解决utf-8不行for file in file_list:if file.endswith('.tif'):print('Processing >>> ' + file)in_raster = gdal.Open(os.path.join(in_dir, file))out_raster = os.path.join(out_dir,file)ds = gdal.Warp(out_raster,# 裁剪图像保存完整路径（包括文件名）in_raster,# 待裁剪的影像format = 'GTiff',# 保存图像的格式cutlineDSName = in_shape, # 矢量文件的完整路径#cutlineWhere=in_shape,# optionally you can filter your cutline (shapefile) based on attribute values#cutlineWhere="FIELD = 'whatever'", # clip specific feature#以shp的范围输出cropToCutline=True, # 将目标图像的范围指定为cutline 矢量图像的范围dstSRS='EPSG:4326', # 参考：WGS84outputBounds=r.bbox,copyMetadata=True,dstNodata = 0)              # select the no data value you likeds=None     #do other stuff with ds object, it is your cropped dataset. in this case we only close the dataset.else:print("This file is not '.tif' file ...")if __name__ == '__main__':timestart = datetime.datetime.now()error_url = []out_path = r'G:\气象数据'shp_path = r'G:\气象数据\裁剪矢量\湖北省.shp'StartTime = '2016'EndTime = '2022'url = 'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/' #月平均数据html=requests.get(url).contentdom = BeautifulSoup(html,"lxml") # 解析html文档#我们只要获取标签a的href值，并加上http://即可构造出下载链接。a_list = dom.findAll("a")  # 找出<a>CHIRPS_lists = [a['href'] for a in a_list if a['href'].endswith('.gz')]  # 找到chirps.gz数据Chirps_Time_Info = []for chirps in CHIRPS_lists:savefile = os.path.join(out_path,chirps)Chirps_Time = re.split(r'[_,.,\s ]\s*', chirps)[-4] #将Chirps文件的时间信息分割if Chirps_Time <= EndTime and Chirps_Time >= StartTime :Chirps_Time_Info.append(chirps)#download(savefile,'https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/'+chirps)IDMdownload('https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_monthly/tifs/'+chirps,savefile,chirps)print("------------------------------------")print("%s文件已添加至下载列表"%chirps)print("------------------------------------")#通过循环判断路径下文件是否下载完成，如果完成则进行Check，若未完成则再次读取路径下文件信息for i in range(1000):Judge = Judge_check(Chirps_Time_Info,out_path)if Judge == True :Check_file(Chirps_Time_Info,out_path)breakelse:time.sleep(120)#120S后重新读取文件夹下文件信息continueungzip(out_path)warpclip(shp_path,os.path.join(out_path,'Chirps'),os.path.join(out_path,'Chirps_clip'))timeend = datetime.datetime.now()print('Running time: %s Seconds'%(timeend-timestart))sys.exit()

这里注意有些库可以直接pip安装,如gzip、shapefile等,但是向gdal库的话建议和python版本匹配,我这里用的python3.9,装的GDAL库版本为3.1.2,如果后续需要对图像进行处理安装shapely、geopandas等库时,要注意有一个版本匹配问题,如果直接pip的话到后面装其他库很有可能会报错…如有需要可以下载我上传的资源，或者私信我发给你~

四、代码结果

最后结果如图所示：第一个图层为代码裁剪得出，第二个图层为arcgis工具（按掩膜提取）得出，第三个为解压缩出的原始数据，可以看到三者像元值一致，因此初步认为代码并无问题。

Python批量下载CHIRPS气象数据并完成解压裁剪等相关推荐

Python脚本使用CDSAPI批量下载ERA5气象数据
使用Python脚本从 Copernicus Climate Data Store (CDS) 检索气象数据具体地,需要检索变量(geopotential.relative_humidity.tem ...
Python 批量下载 ERA-5 Reanalysis 数据
ECMWF大气再分析数据集 ERA-interim 已被 ERA-5 数据集取代,ERA-5 详细信息看这里下面介绍下载的具体步骤: 1. 要下载 ERA-5 数据集,需要先注册一个 CDS 账号. ...
python读取哨兵卫星数据_科学网—【Python】批量下载Sentinel-2卫星数据 - 江佳乐的博文...
简介:本文介绍了一种基于Python批量下载Sentinel-2卫星数据的方法.通过该方法,用户可自定义感兴趣研究区(ROI).限定时间范围.选择产品类别.进行云量筛选,实现自动化批量下载Sentin ...
使用Python批量下载哨兵一号（sentinel-1）的精密轨道数据（precise orbit data）
使用Python批量下载哨兵一号(sentinel-1)的精密轨道数据 1. 安装Python及data_downloader包 2. 下载精密轨道数据与辅助数据 1. 安装Python及data_d ...
用Python批量下载DACC的MODIS数据
本人初次尝试用Python批量下载DACC的MODIS数据,记下步骤,提醒自己,数据还在下载,成功是否未知,等待结果中...... 若有大佬发现步骤有不对之处,望指出,不胜感激. 1.下载Python ...
使用python批量下载天猫数据并进行合并（非爬虫）
使用python批量下载天猫数据并进行合并(非爬虫) 做电商运营少不了数据分析支持,无论是选品.选关键词.研究竞品还是开直通车.店铺引流都需要有数据支持.不过生意参谋虽然强大,但现在都只能显示&quo ...
python下载网页里面所有的图片-Python批量下载网页图片详细教程
很多朋友在网上查找批量下载图片的方法~发觉挺凌乱的,无从下手.这里绿茶小编就来跟大家分享下使用Python批量下载图片方法. 目标:爬取某个网站上n多页的链接,每个链接有n多张图片,每一页对应一个文件 ...
python批量下载b站_python 批量下载bilibili视频的gui程序
运行效果: 完整代码: # !/usr/bin/python # -*- coding:utf-8 -*- # time: 2019/07/02--08:12 __author__ = 'Henry' ...
python怎么批量下载年报_使用Python批量下载Wind数据库中的PDF报告
原标题:使用Python批量下载Wind数据库中的PDF报告背景最近小编出于工作需要,准备在Wind金融数据终端批量下载上市公司2019年第一季度业绩预告.通过相关的条件检索,发现其相关数据有近百 ...

Python批量下载CHIRPS气象数据并完成解压裁剪等

文章目录

前言