Эффективно маскировать и вычислять средние значения для нескольких групп в xarDataset xarray - PullRequest
0 голосов
/ 29 января 2020

У меня есть два xr.Dataset объекта. Один - это непрерывная карта некоторой переменной (здесь precipitation). Другая - это категориальная карта набора регионов ['region_1', 'region_2', 'region_3', 'region_4'].

. Я хочу вычислить среднее значение precip в каждом region в каждом timestep, маскируя по региону / времени и затем выводя кадр данных выглядит следующим образом.

In [6]: df.head()
Out[6]:
    datetime region_name          mean_value
0 2008-01-31    region_1   51.77333333333333
1 2008-02-29    region_1   44.87555555555556
2 2008-03-31    region_1   50.88444444444445
3 2008-04-30    region_1   48.50666666666667
4 2008-05-31    region_1  47.653333333333336

У меня есть код, но он работает очень медленно для реальных наборов данных. Может ли кто-нибудь помочь мне оптимизировать?

Минимальный воспроизводимый пример

Инициализация наших объектов, две переменные одной формы. region объект будет считан из шейп-файла и будет иметь более двух областей.

import xarray as xr
import pandas as pd
import numpy as np

def make_dataset(
    variable_name='precip',
    size=(30, 30),
    start_date='2008-01-01',
    end_date='2010-01-01',
    lonmin=-180.0,
    lonmax=180.0,
    latmin=-55.152,
    latmax=75.024,
):
    # create 2D lat/lon dimension
    lat_len, lon_len = size
    longitudes = np.linspace(lonmin, lonmax, lon_len)
    latitudes = np.linspace(latmin, latmax, lat_len)
    dims = ["lat", "lon"]
    coords = {"lat": latitudes, "lon": longitudes}

    # add time dimension
    times = pd.date_range(start_date, end_date, name="time", freq="M")
    size = (len(times), size[0], size[1])
    dims.insert(0, "time")
    coords["time"] = times

    # create values
    var = np.random.randint(100, size=size)

    return xr.Dataset({variable_name: (dims, var)}, coords=coords), size

ds, size = make_dataset()

# create dummy regions (not contiguous but doesn't matter for this example)
region_ds = xr.ones_like(ds).rename({'precip': 'region'})
array = np.random.choice([0, 1, 2, 3], size=size)
region_ds = region_ds * array

# create a dictionary explaining what the regions area
region_lookup = {
    0: 'region_1',
    1: 'region_2',
    2: 'region_3',
    3: 'region_4',
}

Как выглядят эти объекты?

In[]: ds

Out[]:
<xarray.Dataset>
Dimensions:  (lat: 30, lon: 30, time: 24)
Coordinates:
  * lat      (lat) float64 -55.15 -50.66 -46.17 -41.69 ... 66.05 70.54 75.02
  * lon      (lon) float64 -180.0 -167.6 -155.2 -142.8 ... 155.2 167.6 180.0
  * time     (time) datetime64[ns] 2008-01-31 2008-02-29 ... 2009-12-31
Data variables:
    precip   (time, lat, lon) int64 51 92 14 71 60 20 82 ... 16 33 34 98 23 53

In[]: region_ds

Out[]:
<xarray.Dataset>
Dimensions:  (lat: 30, lon: 30, time: 24)
Coordinates:
  * lat      (lat) float64 -55.15 -50.66 -46.17 -41.69 ... 66.05 70.54 75.02
  * time     (time) datetime64[ns] 2008-01-31 2008-02-29 ... 2009-12-31
  * lon      (lon) float64 -180.0 -167.6 -155.2 -142.8 ... 155.2 167.6 180.0
Data variables:
    region   (time, lat, lon) float64 0.0 0.0 0.0 0.0 0.0 ... 1.0 1.0 1.0 1.0

Текущая реализация

Чтобы рассчитать среднее значение переменной в ds в каждой из областей ['region_1', 'region_2', ...] в region_ds в каждый момент времени, мне нужно l oop по ВРЕМЕНИ и РЕГИОНУ.

I oop по каждому РЕГИОНУ, а затем по каждому TIMESTEP в объекте da. Эта операция довольно медленная, так как набор данных становится больше (больше пикселей и больше временных шагов). Может ли кто-нибудь придумать более эффективную / векторизованную реализацию?

Моя текущая реализация очень медленная для всех регионов и времен, которые мне нужны. Есть ли более эффективное использование numpy / xarray, которое даст мне желаемый результат быстрее?

def drop_nans_and_flatten(dataArray: xr.DataArray) -> np.ndarray:
    """flatten the array and drop nans from that array. Useful for plotting histograms.

    Arguments:
    ---------
    : dataArray (xr.DataArray)
        the DataArray of your value you want to flatten
    """
    # drop NaNs and flatten
    return dataArray.values[~np.isnan(dataArray.values)]

#
da = ds.precip
region_da = region_ds.region
valid_region_ids = [k for k in region_lookup.keys()]

# initialise empty lists
region_names = []
datetimes = []
mean_values = []

for valid_region_id in valid_region_ids:
    for time in da.time.values:
        region_names.append(region_lookup[valid_region_id])
        datetimes.append(time)
        # extract all non-nan values for that time-region
        mean_values.append(
            da.sel(time=time).where(region_da == valid_region_id).mean().values
        )

df = pd.DataFrame(
    {
        "datetime": datetimes,
        "region_name": region_names,
         "mean_value": mean_values,
    }
)

Вывод:

In [6]: df.head()
Out[6]:
    datetime region_name          mean_value
0 2008-01-31    region_1   51.77333333333333
1 2008-02-29    region_1   44.87555555555556
2 2008-03-31    region_1   50.88444444444445
3 2008-04-30    region_1   48.50666666666667
4 2008-05-31    region_1  47.653333333333336

In [7]: df.tail()
Out[7]:
     datetime region_name          mean_value
43 2009-08-31    region_4   50.83111111111111
44 2009-09-30    region_4   48.40888888888889
45 2009-10-31    region_4   51.56148148148148
46 2009-11-30    region_4  48.961481481481485
47 2009-12-31    region_4   48.36296296296296

In [20]: df.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 96 entries, 0 to 95
Data columns (total 3 columns):
datetime       96 non-null datetime64[ns]
region_name    96 non-null object
mean_value     96 non-null object
dtypes: datetime64[ns](1), object(2)
memory usage: 2.4+ KB

In [21]: df.describe()
Out[21]:
                   datetime region_name         mean_value
count                    96          96                 96
unique                   24           4                 96
top     2008-10-31 00:00:00    region_1  48.88984800150122
freq                      4          24                  1
first   2008-01-31 00:00:00         NaN                NaN
last    2009-12-31 00:00:00         NaN                NaN

Любая помощь будет принята с благодарностью! Thankyou

1 Ответ

1 голос
/ 02 февраля 2020

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

regions = xr.concat(
    [(region_ds.region == region_id).expand_dims(region=[region])
     for region_id, region in region_lookup.items()], 
    dim='region'
)
result = ds.precip.where(regions).mean(['lat', 'lon'])

Создает массив данных с размерами 'time' и 'region', где значение в каждой точке является средним значением в данный момент времени в заданном регионе. Было бы просто расширить это значение до средневзвешенного по площади, если бы это тоже было необходимо.


Альтернативный вариант, который дает такой же результат, будет:

regions = xr.DataArray(
    list(region_lookup.keys()),
    coords=[list(region_lookup.values())],
    dims=['region']
)
result = ds.precip.where(regions == region_ds.region).mean(['lat', 'lon'])

Здесь regions - это просто представление DataArray словаря region_lookup.

...