commit
This commit is contained in:
parent
684ae9d4e8
commit
dbcff95cde
|
@ -0,0 +1,36 @@
|
||||||
|
from osgeo import gdal
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
def read_grib_file(file_path):
|
||||||
|
# 打开 grib 文件
|
||||||
|
grib_dataset = gdal.Open(file_path)
|
||||||
|
|
||||||
|
# 获取数据集的基本信息
|
||||||
|
raster_count = grib_dataset.RasterCount
|
||||||
|
x_size = grib_dataset.RasterXSize
|
||||||
|
y_size = grib_dataset.RasterYSize
|
||||||
|
|
||||||
|
# 读取温度数据并计算每个时间步的平均温度
|
||||||
|
temperature_data = []
|
||||||
|
mean_temperature_per_time_step = []
|
||||||
|
for i in range(raster_count):
|
||||||
|
band = grib_dataset.GetRasterBand(i + 1)
|
||||||
|
temperature_data.append(band.ReadAsArray())
|
||||||
|
mean_temperature_per_time_step.append(np.mean(temperature_data[-1]))
|
||||||
|
|
||||||
|
# 关闭数据集
|
||||||
|
grib_dataset = None
|
||||||
|
|
||||||
|
return raster_count, x_size, y_size, temperature_data, mean_temperature_per_time_step
|
||||||
|
|
||||||
|
# 定义文件路径
|
||||||
|
file_path = r"D:\图层\11111.grib"
|
||||||
|
|
||||||
|
# 读取 grib 文件并获取基本信息以及温度数据
|
||||||
|
raster_count, x_size, y_size, temperature_data, mean_temperature_per_time_step = read_grib_file(file_path)
|
||||||
|
|
||||||
|
# 输出结果
|
||||||
|
print("Raster Count:", raster_count)
|
||||||
|
print("X Size:", x_size)
|
||||||
|
print("Y Size:", y_size)
|
||||||
|
print("Mean Temperature per Time Step:", mean_temperature_per_time_step)
|
|
@ -0,0 +1,134 @@
|
||||||
|
from osgeo import gdal
|
||||||
|
# gdal.UseExceptions()
|
||||||
|
# gdal.PushErrorHandler('CPLQuietErrorHandler')
|
||||||
|
import numpy as np
|
||||||
|
from collections import deque
|
||||||
|
import math, os, re
|
||||||
|
|
||||||
|
|
||||||
|
# 日标准:UTC 0:00-24:00;即北京时间 8:00 - 第二天8:00
|
||||||
|
# 降水/潜在蒸发/蒸发:m -> mm :m * 1000 = mm ; sum
|
||||||
|
# 气温:k -> ℃ : k - 273.15 = ℃ ; mean max
|
||||||
|
# 风速: mean
|
||||||
|
|
||||||
|
def get_ori_values(grib, end, start = 1, mode = 't', day = 365):
|
||||||
|
"""
|
||||||
|
读取数据,初步处理
|
||||||
|
"""
|
||||||
|
if grib.RasterCount != day * 24:
|
||||||
|
print("大小错误")
|
||||||
|
return None
|
||||||
|
out = {}
|
||||||
|
for i in range(start, end):
|
||||||
|
band = grib.GetRasterBand(i)
|
||||||
|
if mode == 'k':
|
||||||
|
out[band.GetMetadata()['GRIB_VALID_TIME']] = band.ReadAsArray() * 1000
|
||||||
|
elif mode == "t":
|
||||||
|
out[band.GetMetadata()['GRIB_VALID_TIME']] = band.ReadAsArray() - 273.15
|
||||||
|
elif mode == "s":
|
||||||
|
out[band.GetMetadata()['GRIB_VALID_TIME']] = band.ReadAsArray()
|
||||||
|
# 这个是为了检验数据排序是否出现了问题
|
||||||
|
if not list(out.keys()) == sorted(out):
|
||||||
|
print("顺序错误")
|
||||||
|
return list(out.values())
|
||||||
|
|
||||||
|
"""
|
||||||
|
组合形式
|
||||||
|
降雨量/蒸散发/潜在蒸散发 - k - sum
|
||||||
|
2mu风速/2mv风速 - s - mean
|
||||||
|
2m温度/2m露点温度 - t - mean
|
||||||
|
2m最大温度 - t - max
|
||||||
|
"""
|
||||||
|
|
||||||
|
year = # 年份:2022
|
||||||
|
day = # 当年天数: 365
|
||||||
|
key = # 提取值的时候用的mode 't'温度相关;'k'降雨蒸发相关; 's'风速相关
|
||||||
|
file_name = # 读取的文件名 eg. "降雨量"
|
||||||
|
total_mode = # 统计用的,'mean'平均, 'max'最大, 'sum'求和
|
||||||
|
|
||||||
|
"""目录结构
|
||||||
|
2022
|
||||||
|
2022_降雨量.grib
|
||||||
|
2022_蒸散发.grib
|
||||||
|
...
|
||||||
|
2021
|
||||||
|
2022_降雨量.grib
|
||||||
|
2022_蒸散发.grib
|
||||||
|
...
|
||||||
|
"""
|
||||||
|
|
||||||
|
res = gdal.Open(f".\\ERA5\\原始数据\\{year}\\{year}_{file_name}.grib")
|
||||||
|
out = get_ori_values(res, start = 1, end = res.RasterCount + 1, mode = key, day = day)
|
||||||
|
out = np.stack(out)
|
||||||
|
|
||||||
|
if file_name == "降雨量":
|
||||||
|
"""
|
||||||
|
处理降雨量的时候,额外计算年度暴雨量
|
||||||
|
暴雨定义:24h > 50
|
||||||
|
那么找到所有24h大于50的窗口
|
||||||
|
合并这些窗口后,计算在窗口内的累计值即可
|
||||||
|
"""
|
||||||
|
driver = gdal.GetDriverByName('GTiff')
|
||||||
|
dataset = driver.Create(f'{year}_暴雨量.tif', res.RasterXSize, res.RasterYSize, 1, gdal.GDT_Float64)
|
||||||
|
dataset.SetGeoTransform(res.GetGeoTransform())
|
||||||
|
dataset.SetProjection(res.GetProjection())
|
||||||
|
band = dataset.GetRasterBand(1)
|
||||||
|
|
||||||
|
def calculate_annual_rainfall(rainfall):
|
||||||
|
# Step 1: 找出所有24小时内降雨量大于50mm的时间段
|
||||||
|
windows = []
|
||||||
|
current_window = None
|
||||||
|
rainfall_queue = deque(maxlen=24)
|
||||||
|
rainfall_sum = 0
|
||||||
|
for i, hour_rainfall in enumerate(rainfall):
|
||||||
|
rainfall_queue.append(hour_rainfall)
|
||||||
|
rainfall_sum += hour_rainfall
|
||||||
|
if len(rainfall_queue) == 24:
|
||||||
|
if rainfall_sum > 50:
|
||||||
|
if current_window is None:
|
||||||
|
current_window = [i - 23, i + 1]
|
||||||
|
else:
|
||||||
|
current_window[1] = i + 1
|
||||||
|
else:
|
||||||
|
if current_window is not None:
|
||||||
|
windows.append(current_window)
|
||||||
|
current_window = None
|
||||||
|
rainfall_sum -= rainfall_queue[0]
|
||||||
|
|
||||||
|
if current_window is not None:
|
||||||
|
windows.append(current_window)
|
||||||
|
|
||||||
|
# Step 2: 合并所有重叠的窗口
|
||||||
|
merged_windows = []
|
||||||
|
for start, end in sorted(windows):
|
||||||
|
if merged_windows and start <= merged_windows[-1][1]:
|
||||||
|
merged_windows[-1][1] = max(merged_windows[-1][1], end)
|
||||||
|
else:
|
||||||
|
merged_windows.append([start, end])
|
||||||
|
|
||||||
|
# Step 3: 累计每个合并后的窗口内的降雨量
|
||||||
|
annual_rainfall = 0.0
|
||||||
|
for start, end in merged_windows:
|
||||||
|
annual_rainfall += sum(rainfall[start:end])
|
||||||
|
|
||||||
|
return annual_rainfall
|
||||||
|
|
||||||
|
band.WriteArray(np.apply_along_axis(calculate_annual_rainfall, 0, out))
|
||||||
|
del dataset
|
||||||
|
|
||||||
|
driver = gdal.GetDriverByName('GTiff')
|
||||||
|
dataset = driver.Create(f'{year}_{file_name}.tif', res.RasterXSize, res.RasterYSize, day, gdal.GDT_Float64)
|
||||||
|
dataset.SetGeoTransform(res.GetGeoTransform())
|
||||||
|
dataset.SetProjection(res.GetProjection())
|
||||||
|
|
||||||
|
for i in range(day):
|
||||||
|
band = dataset.GetRasterBand(i + 1)
|
||||||
|
band.SetMetadata({'day': f"{year}{i + 1:0>3d}"})
|
||||||
|
if total_mode == "sum":
|
||||||
|
band.WriteArray(out[i * 24: i * 24 + 24,:,:].sum(axis = 0))
|
||||||
|
elif total_mode == "max":
|
||||||
|
band.WriteArray(out[i * 24: i * 24 + 24,:,:].max(axis = 0))
|
||||||
|
elif total_mode == "mean":
|
||||||
|
band.WriteArray(out[i * 24: i * 24 + 24,:,:].mean(axis = 0))
|
||||||
|
|
||||||
|
del dataset, res, out
|
|
@ -0,0 +1,45 @@
|
||||||
|
from osgeo import gdal
|
||||||
|
import numpy as np
|
||||||
|
|
||||||
|
def process_layer(layer_path):
|
||||||
|
# 读取上传的图层数据
|
||||||
|
layer_data = gdal.Open(layer_path)
|
||||||
|
|
||||||
|
# 处理后的数据可以存储在某个变量中,以便后续计算水源涵养
|
||||||
|
|
||||||
|
# 示例:假设处理后的数据存储在变量 processed_data 中
|
||||||
|
processed_data = ...
|
||||||
|
|
||||||
|
return processed_data
|
||||||
|
|
||||||
|
def calculate_water_conservation(processed_data):
|
||||||
|
|
||||||
|
# 计算后的结果可以存储在某个变量中,例如 result_data
|
||||||
|
|
||||||
|
# 示例:假设计算后的结果存储在变量 result_data 中
|
||||||
|
result_data = ...
|
||||||
|
|
||||||
|
return result_data
|
||||||
|
|
||||||
|
def save_layer(data, output_path):
|
||||||
|
# 将计算结果保存为GeoTIFF格式的文件
|
||||||
|
driver = gdal.GetDriverByName('GTiff')
|
||||||
|
dataset = driver.Create(output_path, data.shape[1], data.shape[0], 1, gdal.GDT_Float64)
|
||||||
|
dataset.SetGeoTransform(layer_data.GetGeoTransform())
|
||||||
|
dataset.SetProjection(layer_data.GetProjection())
|
||||||
|
band = dataset.GetRasterBand(1)
|
||||||
|
band.WriteArray(data)
|
||||||
|
del dataset
|
||||||
|
|
||||||
|
# 上传的图层文件路径
|
||||||
|
uploaded_layer_path = "uploaded_layer.tif"
|
||||||
|
|
||||||
|
# 处理上传的图层数据
|
||||||
|
processed_data = process_layer(uploaded_layer_path)
|
||||||
|
|
||||||
|
# 计算水源涵养
|
||||||
|
water_conservation_data = calculate_water_conservation(processed_data)
|
||||||
|
|
||||||
|
# 保存水源涵养图层
|
||||||
|
output_path = "water_conservation_layer.tif"
|
||||||
|
save_layer(water_conservation_data, output_path)
|
Loading…
Reference in New Issue