编程交流——第2期

目录

主题一：栅格数据的读取、处理和读出

主题二：图形绘制

附件下载

---

主题一：栅格数据的读取、处理和读出

要求：

（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.youkuaiyun.com/download/revised_one/87678338