osgeo：Rasterio

一、基本操作

• 索引：rasterio.open().xy(rows, cols)
• 矩阵维度的转换：reshape_as_raster$\Rightarrow$(维，行，列), reshape_as_image $\Rightarrow$(行，列，维)

import rasterio

# 打开文件
ds = rasterio.open('./GDAL_data/Test_RGB.tif', mode='r')
ds.meta
# ds.driver, ds.name, ds.count, ds.width, ds.height, ds.bounds, mode, dtype, bounds, transform, crs

{'driver': 'GTiff',
 'dtype': 'uint16',
 'nodata': None,
 'width': 2923,
 'height': 2924,
 'count': 7,
 'crs': CRS.from_epsg(32650),
 'transform': Affine(10.0, 0.0, 375870.0,
        0.0, -10.0, 4480780.0)}

# 读取波段数据
ds.read(1), ds.read()

# 空间索引
ds.index(ds.bounds.left + 100, ds.bounds.top - 50), ds.xy(ds.height // 2, ds.width // 2)

((5, 10), (390485.0, 4466155.0))

# 矩阵维度的转换
from rasterio.plot import reshape_as_raster, reshape_as_image
ds = rasterio.open('./GDAL_data/Test_RGB.tif', mode='r').read()
image = reshape_as_image(ds)
ds.shape, image.shape, reshape_as_raster(image).shape

((7, 2924, 2923), (2924, 2923, 7), (7, 2924, 2923))

二、Rasterio作图

• rasterio.plot.show (source,
  with_bounds=True,
  contour=False,
  contour_label_kws=None,
  ax=None,
  title=None,
  transform=None,
  adjust=‘linear’)
  • source: 数组，或栅格影像(r模式) $\Rightarrow$(维， 行， 列)
  • with_bounds: 是否将图像范围更改为图像的空间边界，而不是像素坐标
  • contour: 是否将光栅数据绘制为等高线
  • contour_label_kws: 用于标记轮廓的关键字参数
  • ax: 用于绘图的轴
  • title: 标题
  • transform：Affine, 如果源为数组，则定义仿射变换
  • adjust：如果为RGB图像，则调整的值归一化（0, 1）;如果’线性’，值将根据每个波段的最小/最大值进行调整。如果没有，没有调整将适用
• rasterio.plot.show_hist (source,
  bins=10,
  masked=True,
  title=‘Histogram’,
  ax=None,
  label=None,
  alpha,)
• rasterio.plot.adjust_band (): 调整波段归一化或拉伸

import rasterio as rio
from rasterio.plot import reshape_as_image, reshape_as_raster
import matplotlib.pyplot as plt  # （行，列，维）
ds = rio.open('./GDAL_data/Test_RGBtif')
a = reshape_as_image(ds.read()[:3, :, :])
print(a.shape)
fig, [[ax1, ax2], [ax3, ax4]] = plt.subplots(2, 2, figsize=(12, 10))
ax1.imshow(a)
rio.plot.show(ds, with_bounds=True, ax=ax2, transform=ds.transform)
# rio.plot.show((ds, 1), contour=True, ax=ax2)
rio.plot.show_hist(ds, ax=ax3, label=['1', '2', '3'], alpha=0.5)
rio.plot.show(ds.read()[2:3, :, :], ax=ax4, cmap='cool')
plt.show()

三、地理参考与坐标变换

Rasterio遵循pyproj，并使用dict格式的proj.4语法作为其本机CRS语法

3.1 不同仿射矩阵形式的转换：gdal和rasterio

from rasterio.transform import Affine
ds = rasterio.open('./GDAL_data/Test_RGB.tif', mode='r')
ds2 = gdal.Open('./GDAL_data/Test_RGB.tif')
ds.transform, dsGetGeoTransform(), ds.transform * (1, 1), ds.transform.to_gdal()  # from_gdal

(Affine(10.0, 0.0, 375870.0,
0.0, -10.0, 4480780.0),
(375870.0, 10.0, 0.0, 4480780.0, 0.0, -10.0),
(375880.0, 4480770.0),
(375870.0, 10.0, 0.0, 4480780.0, 0.0, -10.0))

3.2 坐标参考系的不同表示形式：wkt和pyroj

• crs.from_wkt(): ds.GetProjection()
• crs.from_epsg(): 4326
• crs.from_proj4(): A PROJ4 string like “+proj=longlat …”
• crs.to_wkt/to_epsg/to_proj4/to_string/to_authority/to_dict

ds = rio.open('./GDAL_data/Test_RGBtif')
ds2 = gdal.Open('./GDAL_data/Test_RGB.tif')
ds.crs, rasterio.crs.CRS.from_wkt(dsGetProjection())  # ds.crs.wkt= dsGetProjection()

(CRS.from_epsg(32650), CRS.from_epsg(32650))

四、鸟瞰图

import shutil, rasterio
from rasterio.enums import Resampling
path = shutil.copy('./GDAL_data/Test_RGBtif', './GDAL_data/Test_RGB2_quick.tif')
factors = [2, 4, 8, 16]

dsq = rasterio.open(path, 'r+')
dsq.build_overviews(factors, Resampling.average)
dsq.update_tags(ns='rio_overview', resampling='average')
dsq.close()

src = rasterio.open(path, 'r')
[src.overviews(i) for i in src.indexes]
ds = src.read(out_shape=(3, int(src.height / 10), int(src.width / 10)))
ds.shape

plt.imshow(reshape_as_image(ds[:3,:,:]))
plt.show()

五、裁剪影像

import fiona
import rasterio
import rasterio.mask

with fiona.open("tests/data/box.shp", "r") as shapefile:
    shapes = [feature["geometry"] for feature in shapefile]
with rasterio.open("tests/data/RGB.byte.tif") as src:
    out_image, out_transform = rasterio.mask.mask(src, shapes, crop=True)
    out_meta = src.meta
out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("RGB.byte.masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

六、重投影与重采样

• transform_bounds()：将源栅格的边界坐标转换为目标坐标参考系统，沿边缘加密点，以考虑边缘的非线性转换。
• calculate_default_transform()：将边界转换为目标坐标系，如果未提供分辨率，则计算分辨率，并返回目标转换和维度
• reproject(source, destination=None, src_transform=None, gcps=None,
  rpcs=None, src_crs=None, src_nodata=None, dst_transform=None,
  dst_crs=None, dst_nodata=None, dst_resolution=None,
  src_alpha=0, dst_alpha=0,
  resampling=<Resampling.nearest: 0>,
  num_threads=1, init_dest_nodata=True, warp_mem_limit=0,
  **kwargs,)

import numpy as np
import rasterio
from rasterio import Affine as A
from rasterio.warp import reproject, Resampling

with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src:
    src_transform = src.transform

    # Zoom out by a factor of 2 from the center of the source
    # dataset. The destination transform is the product of the
    # source transform, a translation down and to the right, and
    # a scaling.
    dst_transform = src_transform*A.translation(-src.width/0, -src.height/0)*A.scale(0)

    data = src.read()

    kwargs = src.meta
    kwargs['transform'] = dst_transform

    with rasterio.open('/tmp/zoomed-out.tif', 'w', **kwargs) as dst:

        for i, band in enumerate(data, 1):
            dest = np.zeros_like(band)
            reproject(
                band,
                dest,
                src_transform=src_transform,
                src_crs=src.crs,
                dst_transform=dst_transform,
                dst_crs=src.crs,
                resampling=Resampling.nearest)

            dst.write(dest, indexes=i)

# 重采样
from rasterio.enums import Resampling

with rasterio.open("example.tif") as dataset:
    data = dataset.read(
        out_shape=(dataset.height * 2, dataset.width * 2, dataset.count),
        resampling=resampling.bilinear
    )