目录
主题一:栅格数据的读取、处理和读出
要求:
(1)读取附件中的栅格文件,数据内容见附件data.tif;
(2)将栅格中小于等于80的栅格赋值为80,其余值不变;
(3)栅格数据处理完成后,按照原有属性导出。
# 开发人:hwy
from osgeo import gdal
# 1. 读取栅格数据
raster = gdal.Open(r'C:\Users\郝文宇\Desktop\data.tif') # 打开栅格数据
num_bands = raster.RasterCount #获取波段数
print(num_bands)
band = raster.GetRasterBand(1) # 获取第一个波段
arr = band.ReadAsArray() # 将波段读取为 NumPy 数组
# 2. 栅格数据处理
arr[arr <= 80] = 80 # 将小于等于 80 的值赋值为 80
# 3. 导出栅格数据
out_raster = "output.tif" # 设置输出文件名
driver = gdal.GetDriverByName("GTiff") # 获取 GTiff 驱动程序
out_raster_dataset = driver.CreateCopy(out_raster, raster) # 创建输出数据集,属性与输入数据集相同
out_band = out_raster_dataset.GetRasterBand(1) # 获取输出数据集的第一个波段
out_band.WriteArray(arr) # 将处理后的数组写入输出波段
out_band.FlushCache() # 刷新波段缓存
out_band.SetNoDataValue(band.GetNoDataValue()) # 设置输出波段的无数据值
out_raster_dataset.SetProjection(raster.GetProjection()) # 设置输出数据集的投影信息
out_raster_dataset.SetGeoTransform(raster.GetGeoTransform()) # 设置输出数据集的地理变换信息
# -*-coding:utf-8 -*-
"""
作者:XJQ
"""
from osgeo import gdal,gdalconst
import numpy as np
import matplotlib.pyplot as plt
#导入栅格数据并转化为数组(文件路径,数据类型)
def into_dem(image_path,num_format):
image = plt.imread(image_path)
# 将栅格化数据转化为可写矩阵
dataset = gdal.Open(image_path, gdalconst.GA_ReadOnly)
img_width = dataset.RasterXSize
img_height = dataset.RasterYSize
adf_GeoTransform = dataset.GetGeoTransform()
im_Proj = dataset.GetProjection()
img_data = np.array(dataset.ReadAsArray(0, 0, img_width, img_height), dtype=float)
return img_data
#读取参考栅格
def read_img( filename):
dataset = gdal.Open(filename) # 打开文件
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
im_bands = dataset.RasterCount # 波段数
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵,左上角像素的大地坐标和像素分辨率
im_proj = dataset.GetProjection() # 地图投影信息,字符串表示
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
del dataset
return im_geotrans, im_proj
#按照参考样本属性导出栅格
def write_tif(newpath,im_data,im_Geotrans,im_proj,datatype):
if len(im_data.shape)==3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
diver = gdal.GetDriverByName('GTiff')
new_dataset = diver.Create(newpath, im_width, im_height, im_bands, datatype)
new_dataset.SetGeoTransform(im_Geotrans)
new_dataset.SetProjection(im_proj)
if im_bands == 1:
new_dataset.GetRasterBand(1).WriteArray(im_data)
else:
for i in range(im_bands):
new_dataset.GetRasterBand(i+1).WriteArray(im_data[i])
del new_dataset
#对栅格数据下小于80的栅格赋值为80
def tif_do(dem):
for i in range(dem.shape[0]):
for j in range(dem.shape[1]):
if dem[i,j]<=80:
dem[i,j]=80
return
#导入dem数据
dem = into_dem(r'C:\Users\86136\Desktop\附件\附件\data.tif','int')
#进行数据计算
tif_do(dem)
#读取dem属性
im_geotrans,im_proj = read_img(r'C:\Users\86136\Desktop\附件\附件\data.tif')
#导出处理好的栅格数据
write_tif(r'C:\Users\86136\Desktop\附件\附件\1.tif',dem,im_geotrans,im_proj,gdal.GDT_Float32)#输出Q.tif
"""
作者:dwk
"""
import os
import rasterio
# 获取源文件路径和文件名
src_path = r'C:\Users\16370\Desktop\附件\data.tif'
src_dir, src_filename = os.path.split(src_path)
# 构建输出文件路径和文件名
output_path = os.path.join(src_dir, 'output.tif')
# 读取栅格文件
with rasterio.open(src_path) as src:
# 读取栅格数据
data = src.read(1)
# 对数据进行处理
data[data <= 80] = 80
# 导出处理后的栅格数据
with rasterio.open(
output_path,
'w',
driver='GTiff',
width=src.width,
height=src.height,
count=1,
dtype=data.dtype,
crs=src.crs,
transform=src.transform,
) as dst:
dst.write(data, 1)
"""
作者:lk
"""
from osgeo import gdal
import numpy as np
# 读取数据
dataset = gdal.Open('./data.tif')
data = dataset.GetRasterBand(1).ReadAsArray()
# 数据检查
print(np.sum(data < 80))
# 数据处理
data_2 = np.where(data < 80, 80, data)
# 数据检查
print(np.sum(data_2 < 80))
# 数据输出
row, col = data.shape
GeoTransform = dataset.GetGeoTransform() # 获取数据空间参考
GetProjection = dataset.GetProjection()
driver = gdal.GetDriverByName('GTiff')
dst_ds = driver.Create(r'data2.tif', col, row, 1,gdal.GDT_Float32) # 创建数据集
dst_ds.SetGeoTransform(GeoTransform) # 设置空间参考
dst_ds.SetProjection(GetProjection)
dst_ds.GetRasterBand(1).WriteArray(data_2) # 写入数据
# 姓名:ZK
from osgeo import gdal
import numpy as np
# 读取栅格文件
raster = gdal.Open(r'E:\Pythonxx\data.tif')
# 获取栅格信息
geotransform = raster.GetGeoTransform()
projection = raster.GetProjection()
band = raster.GetRasterBand(1)
# 获取栅格数据
data = band.ReadAsArray()
# 将小于等于80的值修改为80
data[data <= 80] = 80
# 导出栅格数据
driver = gdal.GetDriverByName('GTiff')
out_raster = driver.Create('output.tif', raster.RasterXSize, raster.RasterYSize, 1, band.DataType)
out_raster.SetGeoTransform(geotransform)
out_raster.SetProjection(projection)
out_band = out_raster.GetRasterBand(1)
out_band.WriteArray(data)
# 释放资源
out_band.FlushCache()
out_raster.FlushCache()
band.FlushCache()
raster.FlushCache()
主题二:图形绘制
要求:
(1)读取表格数据,为某地降雨数据,数据内容见附件rain.xlsx;
(2)使用python绘制图形,直观表现降雨随时间的变化;
(3)对于图形类型不作要求,美观形象即可,导出图像并上传。
# 开发人:hwy
import pandas as pd
import matplotlib.pyplot as plt
from pylab import mpl
plt.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False
# 读取 Excel 表格数据
data = pd.read_excel(r'C:\Users\郝文宇\Desktop\rain.xlsx', header=0)
# 将“日期”列转换为 datetime 类型,并设置为索引
data['日期'] = pd.to_datetime(data['日期'], format='%Y-%m-%d')
data.set_index('日期', inplace=True)
data['Rainfall'] = pd.to_numeric(data['Rainfall'])
# 绘制折线图
fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(data.index, data['Rainfall'], color='blue', linewidth=1)
print(data['Rainfall'])
# 设置图形属性
ax.set_xlabel('日期', fontsize=14)
ax.set_ylabel('降雨量 (mm)', fontsize=14)
ax.set_title('1995年1月1日到1999年12月5日的逐日降雨量变化图', fontsize=16)
# 设置坐标轴刻度
ax.tick_params(axis='both', which='major', labelsize=12)
# 限制横坐标范围
ax.set_xlim(pd.to_datetime('1995-01-01'), pd.to_datetime('1999-12-05'))
# 保存图形为 PNG 文件
plt.savefig('rain.png', dpi=300, bbox_inches='tight')
# 显示图形
plt.show()
"""
作者:cyf
"""
import pandas as pd
from matplotlib import pyplot as plt
plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']
plt.rcParams["font.sans-serif"] = ["SimHei"]
df = pd.read_excel(r'C:\Users\chai\Desktop\rain.xlsx')
df.dropna(inplace=True)
def date(para):
if type(para) == int:
delta = pd.Timedelta(str(int(para)) + 'days')
time = pd.to_datetime('1899-12-30') + delta
return time
else:
return para
df["日期"] = df["日期"].apply(date)
plt.plot(df["日期"], df["Rainfall"])
plt.xlabel("日期")
plt.ylabel("降水量")
plt.title("降水量趋势图")
plt.tight_layout()
plt.savefig('降水量趋势图', bbox_inches='tight')
plt.show()
附件下载
附件下载地址:https://download.csdn.net/download/revised_one/87678338