python 网格_python 栅格处理利器之Rasterio

本文主要是Automatize data download​automating-gis-processes.github.io

中栅格图像处理的学习笔记，分享给大家，同时也便于自己记忆以及查看，节省时间。数据准备

如果手头有数据，可以跳过该步骤，如果没有，可以按照原文提供的下载链接，保存数据。

import os

import urllib

def get_filename(url):

"""

Parses filename from given url

"""

if url.find('/'):

return url.rsplit('/', 1)[1]

# 保存路径

outdir = r"data"

# 如果没有路径，在当前环境路径下创建

if not os.path.exists(outdir):

os.makedirs(outdir)

# url

url_list = ["https://github.com/Automating-GIS-processes/CSC18/raw/master/data/Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif"]

# 数据下载

for url in url_list:

# Parse filename

fname = get_filename(url)

outfp = os.path.join(outdir, fname)

# Download the file if it does not exist already

if not os.path.exists(outfp):

print("Downloading", fname)

r = urllib.request.urlretrieve(url, outfp)读取

以前面下载的数据为例，数据读取操作如下，rasterio的数据读取还是很python的，这也是分享的主要的原因。

import rasterio

import os

import numpy as np

%matplotlib inline

# 数据路径，根据自己的需求设定

data_dir = "data"

fp = os.path.join(data_dir, "Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif")

#

raster = rasterio.open(fp)

type(raster)

数据读取之后，需要查看栅格的属性，不得不说，感觉这这方面还是Matlab方便，读取之后能直接看到各种属性的变量，直接手动点击就可以方便查看。

地理栅格数据的主要的属性有坐标系、仿射变换(Affine transform)、维度、波段、缺失值格式以及边界等，实现如下

raster.crs # 坐标系

raster.transform # 仿射变换

(raster.width,raster.height) # 维度

raster.count # 波段

raster.nodatavals # 缺失值

raster.dirver # 数据格式

# 上面的所有信息，也可以通过raster.meta一次展示

raster.meta

如果需要查看波段的具体的数值以及统计信息，则可以通过如下的操作实现

band1 = raster.read(1) # 第一波段属性值

# 所有波段

array = raster.read()

# 统计每个波段的信息

stats = []

for band in array:

stats.append({

'min': band.min(),

'mean': band.mean(),

'median': np.median(band),

'max': band.max()})

# 展示统计信息

stats可视化

对于前面的读入的数据，通过以下简单的命令即可

import matplotlib.pyplot as plt

from rasterio.plot import show

%matplotlib inline

# band 1

show((raster, 1))

fig, (ax1, ax2, ax3) = plt.subplots(ncols=3, nrows=1, figsize=(10, 4), sharey=True)

# rgb 展示

show((raster, 4), cmap='Reds', ax=ax1)

show((raster, 3), cmap='Greens', ax=ax2)

show((raster, 1), cmap='Blues', ax=ax3)

# 添加标题

ax1.set_title("Red")

ax2.set_title("Green")

ax3.set_title("Blue")图像合成

假彩色或者真彩色的合成，此处以假彩色为例，真彩色的合成同理。

# 假彩色合成

nir,red,green = raster.read(4),raster.read(3),raster.read(2)

def normalize(array):

"""Normalizes numpy arrays into scale 0.0 - 1.0"""

array_min, array_max = array.min(), array.max()

return ((array - array_min)/(array_max - array_min))

# 标准化

nirn,redn,greenn = normalize(nir),normalize(red),normalize(green)

# 图像合成

nrg = np.dstack((nirn, redn, greenn))

plt.imshow(nrg)

栅格属性的直方图展示

# 栅格

from rasterio.plot import show_hist

show_hist(raster, bins=50, lw=0.0, stacked=False, alpha=0.3,

histtype='stepfilled', title="Histogram")矢量裁剪

栅格中最为常见的功能之一，rasterio中的操作可以通过如下代码实现，直接运行原文的代码可能会因为坐标系的转换出现错误，此处稍微做了修改。

此处涉及到的geopandas包是矢量处理的利器，依赖包包括shapely、fiona等，可以直接通过conda install geopandas 的方式安装，省去麻烦，pycrs包可以通过pip安装。

栅格的保存也是比较重要的操作，在裁剪的过程中保存结果的代码也在此处提及。

from rasterio.mask import mask

from shapely.geometry import box

import geopandas as gpd

from fiona.crs import from_epsg

import pycrs

%matplotlib inline

# 保存

out_tif = os.path.join(data_dir, "Helsinki_Masked.tif")

# 根据raster边界创建边界

minx, miny = 700000,6660000

maxx, maxy = 730000, 6690000

bbox = box(minx, miny, maxx, maxy)

geo = gpd.GeoDataFrame({'geometry': bbox}, index=[0], crs=raster.crs.data)

# geo.plot()

def getFeatures(gdf):

"""Function to parse features from GeoDataFrame in such a manner that rasterio wants them"""

import json

return [json.loads(gdf.to_json())['features'][0]['geometry']]

coords = getFeatures(geo)

# 裁剪

out_img, out_transform = mask(dataset=raster, shapes=coords, crop=True)

# 输出属性信息

out_meta = raster.meta.copy()

out_meta.update({"driver": "GTiff",

"height": out_img.shape[1],

"width": out_img.shape[2],

"transform": out_transform,

"crs":raster.crs}

)

with rasterio.open(out_tif, "w", **out_meta) as dest:

dest.write(out_img)栅格数值运算

以NDVI计算为例，计算的结果的保存可以参照裁剪的计算，并保存ndvi图像。

# red波段

red = raster.read(3)

# 近红外-NIR波段

nir = raster.read(4)

# 数值类型转换

red = red.astype('f4')

nir = nir.astype('f4')

np.seterr(divide='ignore', invalid='ignore')

# ndvi 计算

ndvi = (nir - red) / (nir + red)

#

plt.imshow(ndvi, cmap='terrain_r')

# 添加colorbar

plt.colorbar()

其余的部分在后续进行补充，最近在项目中涉及到的矢量以及栅格的运算的操作，也计划慢慢分享给大家。