Regression with xr.apply_ufunc

[Murk89]

Hi,

I am trying to apply the ordinary least squares model from pythons's sklearn, to every gridbox of my dataset called 'merged'
This error comes up:
ValueError: Expected 2D array, got 1D array instead:
array=[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.].
Reshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.

Would really appreciate some help with this error or an idea of how to apply multivariate regression using xarray's apply_ufunc

from sklearn import linear_model
import salem
import matplotlib
import scipy
import xarray as xr
from salem import mercator_grid, Map, get_demo_file, open_xr_dataset
import pandas as pd
import numpy as np

merged
<xarray.Dataset>
Dimensions:      (year: 20, west_east: 114, south_north: 111)
Coordinates:
  * year         (year) int64 1998 2001 2002 2003 2005 ... 2017 2018 2019 2020
    band         int32 1
  * west_east    (west_east) float64 -2.766e+05 -2.666e+05 ... 8.534e+05
  * south_north  (south_north) float64 -1.353e+05 -1.253e+05 ... 9.647e+05
    lat          (south_north, west_east) float32 62.68 62.68 ... 71.2 71.17
    lon          (south_north, west_east) float32 4.579 4.774 ... 34.04 34.29
Data variables:
    band_data    (year, south_north, west_east) float64 nan nan nan ... nan nan
    SNOWH        (year, south_north, west_east) float32 0.0 0.0 0.0 ... 0.0 0.0
Attributes:
    long_name:     NDVI_mean
    grid_mapping:  spatial_ref
    pyproj_srs:    +proj=lcc +lat_0=64 +lon_0=10 +lat_1=64 +lat_2=68 +x_0=0 +...


def regress (x,y):
 ols= linear_model.LinearRegression(),
 model.score,model.coef_ ,model.intercept_ = ols.fit(x, y),
 return xr.dataset([model.score,model.coef_ ,model.intercept_])

trial = xr.apply_ufunc(regress, merged['SNOWH'], merged['band_data'],
                       input_core_dims=[['year'], ['year']],
                       vectorize=True,
                       dask="parallelized",
                       )

[max-sixty]

Hi @Murk89 — could you rework the example consistent with #5404? Thank you

[alessioarena]

Hi @Murk89
There is nothing wrong on the xarray side. vectorize=True will basically iterate over the non input_core_dims and pass them to your function as 1D array.
However, from sklearn documentation seems that you need to pass "{array-like, sparse matrix} of shape (n_samples, n_features)", so in your case I suspect that as a minimum you need to provide an array of shape (1, 20) (2D array with 1 sample and 20 features).

To do that modify this part of the code:

def regress (x,y):
 ols= linear_model.LinearRegression(),
 model = ols.fit(x.reshape(1, -1), y.reshape(1, -1)),
 return [model.score, model.coef_ , model.intercept_]

Hope this fixes your issue

[Murk89]

Hi @alessioarena,
Thanks for your help. I tried with your suggestion. Now get this error:
ValueError: not enough values to unpack (expected 3, got 1)
Any idea about this?

[alessioarena]

Updated my answer. fit returns the model, than you can access the model attributes separately, like intercept_.

Quick note, you will also have to change some parameters in apply_ufunc. Basically you need to explicitly tell to apply_ufunc how the result will look like. You can do that by setting:

• output_core_dims, to set the name of the new axis, (e.g. output_core_dims = [['model_characteristics']])
• output_dtypes, to set the data type (e.g. output_dtypes = ['float64'])
• output_sizes, to set the size of the new dimension (e.g. dask_gufunc_kwargs={'output_sizes': {'model_characteristics': 3}})
  You may have to fiddle around with few other parameters, but if you follow the error messages and this guide you should be able to set the rest

[Illviljan]

Wouldn't polyfit or curve_fit work for this?