“栅格”还起源于电视技术，一种扫描模式（如阴极射线管中的电子束），其中从上到下成行从一侧到另一侧扫描一个区域。
从栅格这个词的来源，我们可以看出耕地和电视扫描，都是在进行一种网格化的过程。
这种网格化的结果：
产生了由像素（微小的彩色方块）网格构建的图像，这种图像就是栅格图像。
Google 地球等工具中查看过照片或图像，那么您之前就已经查看过并使用过栅格。
但是，栅格数据与照片不同，因为它们是空间参考的。每个像素代表地面上的一块土地。
地理空间栅格与数码照片的唯一不同之处在于它将数据连接到特定位置的空间信息。
创建两个ndarray对象，一个X跨越 [-90°,90°] 经度，并Y覆盖 [-90°,90°] 纬度。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('TkAgg')
x = np.linspace(-90, 90, 6)
y = np.linspace(90, -90, 6)
X, Y = np.meshgrid(x, y)
plt.figure(figsize=(6, 4))
# 散点图表示网格坐标
plt.scatter(X, Y)
# 绘制垂直线
for i in range(len(x)):
plt.vlines(x[i], ymin=min(Y.flatten()), ymax=max(Y.flatten()), colors='r', linestyles='-.', linewidth=0.5)
# 绘制水平线
for i in range(len(y)):
plt.hlines(y[i], xmin=min(X.flatten()), xmax=max(X.flatten()), colors='g', linestyles='-.', linewidth=0.5)
plt.show()
View Code
生成温度栅格数据：
import matplotlib.pyplot as plt
import matplotlib
import numpy as np
matplotlib.use('TkAgg')
x = np.linspace(-90, 90, 6)
y = np.linspace(90, -90, 6)
X, Y = np.meshgrid(x, y)
Z1 = np.abs(((X- 10) ** 2 + (Y - 10) ** 2) / 1 ** 2)
Z2 = np.abs(((X+ 10) ** 2 + (Y + 10) ** 2) / 2.5 ** 2)
Z = (Z1 - Z2)
plt.imshow(Z)
plt.title("Temperature")
for i in range(len(x)):
plt.vlines(x[i], ymin=min(Y.flatten()), ymax=max(Y.flatten()), colors='r', linestyles='-.', linewidth=0.5)
# 绘制水平线
for i in range(len(y)):
plt.hlines(y[i], xmin=min(X.flatten()), xmax=max(X.flatten()), colors='g', linestyles='-.', linewidth=0.5)
plt.show()
View Code
生成了z数组，它具有6行、6列。
当我们不对这组数据做任何处理，可以看出，这组数据处于网格的中心位置，Z不包含有关其位置的数据。
它目前的位置是它本身的行列号，比如放大上图，可以看到，它的左上角位置为（0，0）。
这与在python里查看数组z，它的行列号索引一样。
有了数组z，这是数据。
生成了一个网格，这个网络的一些坐标存储在（X，Y）。
需要将数据与网格建立联系：
- 定义左上角的坐标，比如我们可以定义数据的左上角坐标为（90°，90°），就像下面程序：Affine.translation(x[0], y[0])
- 定义单元分辨率，也就是数据的1格和网络上的1格之间的比例关系，比如我们定义数据的1格就是1°，那么我们就可以继续写程序，Affine.translation(x[0], y[0]) * Affine.scale(1, 1)
- 实现数据和网络的关联，通过rasterio库创建一个栅格数据，把它的transform设置为我们定义好的左上角坐标和单元分辨率。
rio.open()函数打开或创建一个 Rasterio 数据集。在这里，它创建一个新的 GeoTIFF 文件用于写入数据。driver="GTiff"：指定输出文件类型为 GeoTIFF。height=Z.shape[0]和width=Z.shape[1]：指定输出栅格数据的高度和宽度，即数组 Z 的形状。count=1：指定数据集中的波段数量，这里是一个单波段数据。dtype=Z.dtype：指定输出数据的数据类型，与数组 Z 的数据类型相同。crs="+proj=latlong"：指定输出数据的坐标参考系统（Coordinate Reference System，CRS），这里是经纬度投影。transform=transform1：定义输出栅格数据的空间参考，包括位置和分辨率。
from rasterio.transform import Affine import rasterio as rio import numpy as np x = np.linspace(-90, 90, 6) y = np.linspace(90, -90, 6) X, Y = np.meshgrid(x, y) Z1 =np.abs(((X - 10) ** 2 + (Y - 10) ** 2) / 1 ** 2) Z2 =np.abs(((X + 10) ** 2 + (Y + 10) ** 2) / 2.5 ** 2) Z =(Z1 - Z2) res = (x[-1] - x[0]) / 240.0 res=10 transform = Affine.translation(x[0] - res / 2, y[0] - res / 2) * Affine.scale(res, res) transform1 = Affine.translation(x[0], y[0]) * Affine.scale(1, 1) # open in 'write' mode, unpack profile info to dst from pathlib import Path folder_path = Path('D:/modis/') file_name = 'new_raster.tif' full_path = folder_path / file_name # 检查路径是否存在 if not folder_path.exists(): folder_path.mkdir(parents=True) # 如果文件夹不存在，创建文件夹 with rio.open( "D:/python/modis/new_raster3.tif", "w", driver="GTiff", # output file type height=Z.shape[0], # shape of array width=Z.shape[1], count=1, # number of bands dtype=Z.dtype, # output datatype crs="+proj=latlong", # CRS transform=transform1, # location and resolution of upper left cell ) as dst: # check for number of bands if dst.count == 1: # write single band dst.write(Z, 1) else: # write each band individually for band in range(len(Z)): # write data, band # (starting from 1) dst.write(Z[band], band + 1)View Code