После разговора о чтении и сохранении данных в формате Geotiff в этой статье речь пойдет о том, как использовать Python для обработки серии растровых данных (в этой статье в качестве примера приведены временные ряды).

Предположим, у нас есть ряд осадков каждого года в определенной области.,Десятилетия в общей сложности,хочу этого сейчасИзменение тренда годовых осадков в каждом пикселе и проверка значимости тренда (получение значения P),Как это сделать?

Идеи

Для растровых данных,Он включает в себяМетаинформация + данные。Находим тренд изменения годового количества осадков для каждого пикселя и соответствующийPценить,На самом деле просто обработка данных,Метаинформация практически не изменяется.

в процессе обработки,Мы хотим, чтобы каждый пиксель находился внутриизмерение временименяющиеся тенденции в,Аналогично изображению ниже:

Цитата с сайта ArcGIS

То есть для каждой из приведенных выше цифр мыВременной ряд бараПросто ищите тенденции。Понятно Идеи,Это очень просто,Перейдем непосредственно к коду.

Создание данных

Создание сказано здесь данных То есть объединить наши несколько растровых последовательностей во что-то похожее на приведенное выше.Куб пространства-времени (здесь относится только к сетке пространства-времени)。

Ранее мы говорили о том, как читатьодиночная сетка,После прочтения это пустой ndarray,Затем просто выполните соответствующее сращивание матриц:

1. Импортировать пакет

import rasterio
import scipy.stats  as ss 
import numpy as np
import glob
from rasterio.plot import show

1. Прочтите метаданные первой картинки, чтобы облегчить запись результата в конце.

with rasterio.open(rs_files[0]) as src:
    show(src)
    profile = src.profile
    data = src.read()

1. сформировать куб пространства-времени

Здесь также очень просто сформировать растровый куб. Просто используйте функцию конкатенации numpy, чтобы объединить его:

ds = [rasterio.open(f).read() for f in rs_files]
ds = np.concatenate(ds,axis=0)
ds.shape

>> (36, 133, 110)

show(ds[1,:,:])

Расчет тренда и значения p

Как упоминалось ранее, пока каждый Временной ряд Просто используйте бара для расчета тренда.,Так как же посчитать все Временной ряд бара одновременно?Великие боги приготовили для нас соответствующие инструменты.,Это функция apply_along_axis numpy.,Подробности смотрите по ссылке [2].

Проще говоря, эта функция можетПримените функцию, которую мы определили, вдоль определенного измерения.。Поэтому нам нужно только определить соответствующую функцию,вдольизмерение Применить:

1. определить функцию

def rs_slope(y,method='linear'):
    """
    Compute the slope of a series of raster data along a given axis(always be time).
    """
    l = len(y)
    ys = range(l)
    try:
        if method == 'linear':
            slope, intercept, r_value, p_value, std_err = ss.linregress(ys,y)
    except:
        slope = np.nan
    return slope, r_value**2, p_value

1. функция приложения

slope_rs,r2,pv = np.apply_along_axis(rs_slope, 0,arr=ds)

1. Сохранить результаты

profile['nodata'] = 0
with rasterio.open('slope.tif','w',**profile) as dest:
    dest.write(slope_rs,1)

1. Посмотреть результаты

with rasterio.open('./slope.tif') as src:
    show(src)

На этом этапе расчет линейного тренда для каждого пикселя завершен. Однако приведенный выше код сохраняет только значение тренда и не сохраняет квадрат R и значение p. Читатель может изменить его в соответствии с кодом.

небольшая добавка

Почему он не писался год? Потому что при расчете тренда,Если вас не волнует перехват,Так год 0-35 или 1980-2015?,Рассчитанное вами значение тренда (т. е. a в формуле ниже),х — год) одинаковы,Тогда нет необходимости тратить больше вычислительной мощности:

Подвести итог

1. При обработке растровых последовательностей метаинформация, как правило, остается неизменной, поэтому метаинформацию определенных исходных данных можно использовать в качестве шаблона для облегчения сохранения обработанных результатов;
2. Что касается значения растровых данных, это просто массив. Умело используя функции numpy, мы можем достичь многих необходимых нам функций;
3. Если алгоритм можно упростить, упростите его настолько, насколько это возможно.

ссылка

【1】

https://pro.arcgis.com/zh-cn/pro-app/latest/tool-reference/big-data-analytics/create-space-time-cube.htm

【2】https://numpy.org/doc/stable/reference/generated/numpy.apply_along_axis.html