Hướng dẫn raster calculation in python - tính toán raster trong python

Tiến hành tính toán giữa các dải hoặc raster là một nhiệm vụ GIS phổ biến khác. Ở đây, chúng tôi sẽ tính toán NDVI (Chỉ số thực vật khác biệt được chuẩn hóa) dựa trên bộ dữ liệu Landsat mà chúng tôi đã tải xuống từ Vùng Helsinki. Tiến hành tính toán với rasterio khá đơn giản nếu mức độ, vv khớp với các giá trị của các raster được lưu trữ dưới dạng mảng numpy (tương tự như các cột được lưu trữ trong Geo/pandas, tức là ____10).

Hãy bắt đầu bằng cách nhập các mô -đun cần thiết

# Filepath
In [4]: fp = r"C:\HY-DATA\HENTENKA\CSC\Data\Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif"

In [5]: raster = rasterio.open(fp)

1 và numpy

In [1]: import rasterio

In [2]: import numpy as np

In [3]: from rasterio.plot import show

Hãy để đọc các tập tin mà chúng tôi đeo mặt nạ cho vùng Helsinki

# Filepath
In [4]: fp = r"C:\HY-DATA\HENTENKA\CSC\Data\Helsinki_masked_p188r018_7t20020529_z34__LV-FIN.tif"

In [5]: raster = rasterio.open(fp)

Để tính toán NDVI (Chỉ số thực vật khác biệt được chuẩn hóa), bạn cần hai dải: Band-4 là Kênh Đỏ và Band-5 là hồng ngoại gần (NIR)

Hãy để đọc các ban nhạc từ nguồn raster của chúng tôi

In [6]: red = raster.read(4)

In [7]: nir = raster.read(5)

In [8]: red
Out[8]: 
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

In [9]: nir
Out[9]: 
array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

In [10]: type(nir)
Out[10]: numpy.ndarray

In [11]: show(nir)
Out[11]: <matplotlib.axes._subplots.AxesSubplot at 0x20c905720f0>

Như chúng ta có thể thấy các giá trị được lưu trữ dưới dạng numpy.ndarray.

Hãy để thay đổi loại dữ liệu từ UINT8 sang nổi để chúng ta có thể có các số điểm nổi được lưu trữ trong các mảng của chúng ta.

In [12]: red = red.astype(float)

In [13]: nir = nir.astype(float)

In [14]: nir
Out[14]: 
array([[ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       ..., 
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

Bây giờ chúng ta có thể thấy rằng các số đã thay đổi thành số thập phân (có một dấu chấm sau 0).

Tiếp theo chúng ta cần điều chỉnh hành vi của Numpy một chút. Theo mặc định, Numpy sẽ phàn nàn về việc chia với các giá trị bằng không. Chúng ta cần thay đổi hành vi đó vì chúng ta có nhiều giá trị 0 trong dữ liệu của mình.

Cho phép 0 phân chia trong Numpy

In [15]: np.seterr(divide='ignore', invalid='ignore')
Out[15]: {'divide': 'warn', 'invalid': 'warn', 'over': 'warn', 'under': 'ignore'}

Bây giờ chúng ta cần khởi tạo NDVI với số không trước khi chúng ta thực hiện các tính toán (đây là thủ thuật cụ thể không có

In [16]: ndvi = np.empty(raster.shape, dtype=rasterio.float32)

Bây giờ chúng tôi đã sẵn sàng để tính toán NDVI. Đầu tiên, chúng ta có thể tạo một bộ lọc trong đó chúng ta tính toán các giá trị trên các pixel đó có giá trị lớn hơn 0.

In [17]: check = np.logical_or ( red > 0, nir > 0 )

Bây giờ chúng ta có thể áp dụng bộ lọc đó và tính toán chỉ số NDVI.

In [18]: ndvi = np.where ( check,  (nir - red ) / ( nir + red ), -999 )

In [19]: ndvi
Out[19]: 
array([[-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       ..., 
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.],
       [-999., -999., -999., ..., -999., -999., -999.]])

In [20]: ndvi.mean()
Out[20]: -108.16964389074909

In [21]: ndvi.std()
Out[21]: 310.4604815144337

In [22]: show(ndvi, cmap='summer')
Out[22]: <matplotlib.axes._subplots.AxesSubplot at 0x20c99ecd780>