目录
- 结果
- 输入文件
- 核心代码
结果
输入文件
1、需要有经纬度信息以及对应的单点值
2、还要用到一个研究区的矢量文件,当然上面点的经纬度信息要在该矢量文件以内
核心代码
file_path = workspace1
# Attempt to read the Excel file
df = readDataFile(file_path)
date = str("Value")
df_cleaned = df.dropna(subset=[date])
shp_file_path = workspace2
region_shp = gpd.read_file(shp_file_path)
bounds = region_shp.bounds
data = pd.read_excel(file_path)
coordinates = data[['POINT_X', 'POINT_Y']].values
values = data[date].values
longitude = data['POINT_X'].values
latitude = data['POINT_Y'].values
OK = OrdinaryKriging(longitude, latitude, values, variogram_model='spherical', verbose=False, enable_plotting=False)
grid_x, grid_y = np.mgrid[bounds.minx.min():bounds.maxx.max():400j, bounds.miny.min():bounds.maxy.max():400j]
grid_z, ss = OK.execute('grid', grid_x[:, 0], grid_y[0, :])
transform = rasterio.transform.from_origin(grid_x.min(), grid_y.max(), (grid_x.max() - grid_x.min()) / 400,
(grid_y.max() - grid_y.min()) / 400)
out_shape = (grid_z.shape[0], grid_z.shape[1])
mask = geometry_mask(region_shp.geometry, invert=True, transform=transform, out_shape=out_shape)
grid_z_masked = np.where(mask, grid_z, np.nan)
transform = from_origin(grid_x.min(), grid_y.max(), (grid_x.max() - grid_x.min()) / 400,
(grid_y.max() - grid_y.min()) / 400)
with rasterio.open(workspace3, 'w', driver='GTiff',
height=grid_z.shape[0], width=grid_z.shape[1],
count=1, dtype=str(grid_z.dtype),
crs='+proj=latlong', transform=transform) as dst:
dst.write(grid_z_masked, 1)
del dst
filename = workspace3
with rasterio.open(filename) as src:
# 读取第一个波段的数据
band_data1 = src.read(1)
# 归一化处理
_normalized = np.floor(255 * (band_data1 - np.nanmin(band_data1)) / (np.nanmax(band_data1) - np.nanmin(band_data1)))
# 应用色彩查找表(LUT)
lut = np.zeros((256, 3), dtype=np.uint8)
# [填充LUT,与原代码相同的逻辑]
for i in range(256):
if i < 51:
# 红色到黄色(增加绿色)
lut[i, 0] = i * 5
lut[i, 1] = 150+i * 2
elif i < 102:
# 黄色到绿色(减少红色)
lut[i, 0] = 255 - (i - 51) * 5
lut[i, 1] = 255
elif i < 153:
# 绿色到青色(增加蓝色)
lut[i, 1] = 255
lut[i, 2] = (i - 102) * 5
elif i < 204:
# 青色到蓝色(减少绿色)
lut[i, 1] = 255 - (i - 153) * 5
lut[i, 2] = 255
else:
# 蓝色到紫色(增加红色)
lut[i, 0] = (i - 204) * 5
lut[i, 2] = 255
# 应用LUT生成伪彩色图像
colored_image = lut[np.int16(_normalized)]
for i in range(3):
colored_image[:, :, i] = np.where(_normalized > 0, colored_image[:, :, i], 0)
# 创建一个新的TIFF文件来保存伪彩色图像
with rasterio.open(
workspace4,
'w',
driver='GTiff',
height=colored_image.shape[0],
width=colored_image.shape[1],
count=3,
dtype=colored_image.dtype,
crs=src.crs,
transform=src.transform,
) as dst:
for i in range(3):
dst.write(colored_image[:, :, i], i + 1)