How to do Spatial correlation analysis and significance value plotting as small dot over the same spatial plot

[HARSHBARDHANK94]

import xarray as xr
import glob
import numpy as np
from scipy import stats

# Fetching .nc files from folders................................................
filename1= glob.glob("E:\\ocean_data\\mld\\*.nc")
filename2= glob.glob("E:\\ocean_data\\sst\\*.nc")

# Reading multiple files using xarray........................................
data1=xr.open_mfdataset(filename1)
data2=xr.open_mfdataset(filename2)

# Area of Interest sub-setting....................................................
AOI1=data1.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))
AOI2=data2.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))

# Data value extraction..................................................................
MLD=AOI1.mixed_layer_depth.values   # shape of (18,26,62)
SST=AOI2.surface_temperature.values # shape of (18,26,62)
lat=AOI1.lat.values  	            # shape of (26,)
lon=AOI1.lon.values                 # shape of (62,)
time=AOI1.time.values               # shape of (18,)

# Linear Reshaping along time.......................................................  
value_1=MLD.reshape(len(time),-1)   # shape of (18,1736)
value_2=SST.reshape(len(time),-1)   # shape of (18,1736)

# Empty List and Array..................................................................
array1= np.array([])
list1=[]
array2= np.array([])
list2=[]

# Looping through the values and apply stats.pearsonr............................................................
for i in range(0,1736):
		ds_1=np.where(np.isfinite(value_1[:,i]),value_1[:,i], np.nanmean(value_1[:,i]))
		ds_2=np.where(np.isfinite(value_2[:,i]),value_2[:,i], np.nanmean(value_2[:,i]))
		corr, p_value=stats.pearsonr(ds_1,ds_2)
		list1.append(corr), list2.append(p_value)
		
array1=np.append(array1, list1)
array1=array1.reshape(28,62)     # shape of (28,62)

array2=np.append(array2, list2)
array2=array2.reshape(28,62)    # shape of (28,62)

# Plotting the array1................................................................................................

''' Note here that array1 contains the correlation value between two data , while
    array2 contains the respective significance (p_value).
'''

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
fig= plt.figure(figsize=(24,12))
ax1 = fig.add_subplot(1,2,1,  projection=ccrs.PlateCarree())
ax2 = fig.add_subplot(1,2,2,  projection=ccrs.PlateCarree())

ax1.pcolormesh(lon,lat, array1 ,cmap=plt.cm.get_cmap('rainbow'),
			   transform=ccrs.PlateCarree())

ax2.pcolormesh(lon,lat, array2 ,cmap=plt.cm.get_cmap('rainbow'),
			   transform=ccrs.PlateCarree())

[HARSHBARDHANK94]

The above shown script , I have written for correlation analysis between two spatial datasets.
Q1.>Does this an appropriate way of doing spatial data correlation ? or some other/way in xarray direct function are available to do it

Q2.> The above script produce two separate plots for r_value and p_value, but I am looking to produce single spatial plot having r_value as plt.pcolormesh plot and p_value<0.5 as small dot over it, wherever applicable !!

Any help from the community !! Thanks in anticipation !!

[HARSHBARDHANK94]

Example plot: (Note: This plot has been taken from a literature and has been produced using MATLAB, as mentioned there!)

[ahuang11]

Maybe relevant correlations:
https://xarray.pydata.org/en/stable/generated/xarray.corr.html
https://xarray.pydata.org/en/stable/user-guide/dask.html?highlight=spearman#apply-ufunc
https://xskillscore.readthedocs.io/en/stable/quick-start.html#Correlation-Based

One way is to first mask with ds.where(p_value < 0.05) and then plot scatter() and scale the size based off the value like this:
https://matplotlib.org/stable/gallery/lines_bars_and_markers/scatter_masked.html#sphx-glr-gallery-lines-bars-and-markers-scatter-masked-py

Another way potential way: https://stackoverflow.com/questions/44250350/how-to-add-stippling-or-hatching-to-individual-polygons-in-a-basemap

[ahuang11]

** This doesn't include the p_value; just showing how to plot it. Use xskillscore if you want p-value.

import xarray as xr
import hvplot.xarray
import xarray as xr

da = xr.tutorial.open_dataset('air_temperature')["air"].load()
da_2013 = da.sel(time="2013")
da_2014 = da.sel(time="2014")
da_2013["time"] = da_2013["time"].dt.strftime("%m%d%H")
da_2014["time"] = da_2014["time"].dt.strftime("%m%d%H")
da_corr = xr.corr(da_2013, da_2014, dim="time")

da_corr.hvplot("lon", "lat", coastline=True, cmap="Reds") * (da_corr.where(da_corr > 0.8)).hvplot.points(
    "lon", "lat", c="air", alpha="air", coastline=True, cmap="Gray").opts("Points", clipping_colors={'NaN': 'transparent'})

[HARSHBARDHANK94]

** This doesn't include the p_value; just showing how to plot it. Use xskillscore if you want p-value.

import xarray as xr
import hvplot.xarray
import xarray as xr

da = xr.tutorial.open_dataset('air_temperature')["air"].load()
da_2013 = da.sel(time="2013")
da_2014 = da.sel(time="2014")
da_2013["time"] = da_2013["time"].dt.strftime("%m%d%H")
da_2014["time"] = da_2014["time"].dt.strftime("%m%d%H")
da_corr = xr.corr(da_2013, da_2014, dim="time")

da_corr.hvplot("lon", "lat", coastline=True, cmap="Reds") * (da_corr.where(da_corr > 0.8)).hvplot.points(
    "lon", "lat", c="air", alpha="air", coastline=True, cmap="Gray").opts("Points", clipping_colors={'NaN': 'transparent'})

Didn't worked well, as expected.
and scatter plot will not work here, bcz it takes 'x' and 'y' positional value having same size, while in our case positional value are lat and lon having different size.
After use I thing hvplot is not accustomed to be used with climatological multiple data plots , that usually requires sub-plotting, formatting and adjustment. But in such cases, hvplot seems numb with its own style, it also takes too much time.

I just wanted to put a simple dot, where p_value<0.05 on trend plot, like:

[HARSHBARDHANK94]

As suggested above, I am using xskillscore to compute the Pearson correlation between two same dimension datasets as earlier
(time: 216, lat:28, lon:62). The Example code below:

import glob
import numpy as np
import xarray as xr
import xskillscore as xs
import matplotlib.pyplot as plt
import dask

# Fetching .nc files from folders...................
filename1= glob.glob("E:\\ACI_data\\AOD\\*.nc")
filename2= glob.glob("E:\\ACI_data\\CER\\*.nc")

# Reading multiple files using xarray
data1=xr.open_mfdataset(filename1)
data2=xr.open_mfdataset(filename2)

# Area of Interest sub-setting.............................................
AOI1=data1.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))
AOI2=data2.sel(lat=slice(-1,26.5), lon=slice(39.5,100.5))

# Selecting variables from the dataarray...................................
AOD=AOI1.MYD08_M3_6_1_AOD_550_Dark_Target_Deep_Blue_Combined_Mean_Mean
CER=AOI2.MYD08_M3_6_1_Cloud_Effective_Radius_Liquid_Mean_Mean

r = xs.pearson_r(AOD, CER, dim='time', skipna=True)

[HARSHBARDHANK94]

But I am getting Value Error:

Traceback (most recent call last):

  File "C:\Users\HARSHB~1\AppData\Local\Temp/ipykernel_1668/4120287431.py", line 1, in <module>
    r=xs.pearson_r(AOD, CER, dim='time', skipna=True)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xskillscore\core\deterministic.py", line 206, in pearson_r
    return xr.apply_ufunc(

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\computation.py", line 1174, in apply_ufunc
    return apply_dataarray_vfunc(

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\computation.py", line 293, in apply_dataarray_vfunc
    result_var = func(*data_vars)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\computation.py", line 681, in apply_variable_ufunc
    raise ValueError(


ValueError: dimension time on 0th function argument to apply_ufunc with 
dask='parallelized' consists of multiple chunks, but is also a core dimension. 
To fix, either rechunk into a single dask array chunk along this dimension, i.e., 
``.chunk(dict(time=-1))``, or pass ``allow_rechunk=True`` in ``dask_gufunc_kwargs`` but
 beware that this may significantly increase memory usage.

Any help on how to set dask_gufunc_kwargs allow_rechunk=True, in the above code.
I am beginner to xarray.apply_ufunc and dask.array.apply_gufunc, as I got these words from previous discussion.

[gkb999]

Are you successful at this stage?
I frequently get memory issues and I'm not successful either at using dask chunks.

[HARSHBARDHANK94]

I tried to write something like that:

def correlation_gufunc(x,y):
	return xs.pearson_r(x, y) 

def pearson_correlation(x, y, dim):
    return xr.apply_ufunc(
        correlation_gufunc,
        x,
        y,
        input_core_dims=[[dim], [dim]],
        dask="parallelized",
        output_dtypes=[float],
		dask_gufunc_kwargs={'allow_rechunk': True},
	)
    
r=pearson_correlation(AOD, CER, "time").compute()

But again I got an another type of error:

Traceback (most recent call last):

  File "C:\Users\HARSHB~1\AppData\Local\Temp/ipykernel_1668/2321186493.py", line 15, in <module>
    r=pearson_correlation(AOD, CER, 'time').compute()

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\dataarray.py", line 955, in compute
    return new.load(**kwargs)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\dataarray.py", line 929, in load
    ds = self._to_temp_dataset().load(**kwargs)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xarray\core\dataset.py", line 865, in load
    evaluated_data = da.compute(*lazy_data.values(), **kwargs)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\base.py", line 570, in compute
    results = schedule(dsk, keys, **kwargs)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\threaded.py", line 79, in get
    results = get_async(

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\local.py", line 507, in get_async
    raise_exception(exc, tb)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\local.py", line 315, in reraise
    raise exc

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\local.py", line 220, in execute_task
    result = _execute_task(task, data)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\core.py", line 119, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\core.py", line 119, in <genexpr>
    return func(*(_execute_task(a, cache) for a in args))

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\core.py", line 119, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\optimization.py", line 969, in __call__
    return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\core.py", line 149, in get
    result = _execute_task(task, cache)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\dask\core.py", line 119, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))

  File "C:\Users\HARSHB~1\AppData\Local\Temp/ipykernel_1668/2321186493.py", line 2, in correlation_gufunc
    return xs.pearson_r(x, y)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xskillscore\core\deterministic.py", line 199, in pearson_r
    dim, _ = _preprocess_dims(dim, a)

  File "C:\Users\HARSHBARDHAN\anaconda3\lib\site-packages\xskillscore\core\utils.py", line 29, in _preprocess_dims
    dim = list(a.dims)

AttributeError: 'numpy.ndarray' object has no attribute 'dims'

[dcherian]

This error means that function xs.pearson_r expects DataArrays not numpy arrays, so it cannot be used with apply_ufunc.
Try rechunking manually

r = xs.pearson_r(AOD.chunk({"time": -1}), CER.chunk({"time": -1}), dim='time', skipna=True)

[aaronspring]

same answer here linking to xarray-contrib/xskillscore#357

[HARSHBARDHANK94]

Thanks a lot, finally:....

[HARSHBARDHANK94]

Still I am looking some robust way to include significance plot as dot mark (p<0.05), wherever applicable over spatial trend plot. Suggestion regarding hvplot.points, seems not compatible with the grouped sub-plots (e.g., 3row*3col=9 plots) . hvplot can be used with single plot. It also takes enough memory and time to render the single plot.

[Asheef1447]

Hi there, were you able to find a more robust way of doing this?

[fllamorim]

Hi, I hope this helps. It is a lot of brute force, but it solved my problem.

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

Correlation

anomalies_chl ('time: 276,lat: 288,lon: 648)

anomalies_sst ('time: 276,lat: 288,lon: 648)

result_cor=[]
p_value=[]

for lat in range(0,len(chl_anomalies.lat)):
for lon in range(0,len(chl_anomalies.lon)):
try:
slope, intercept, r_value, pp, std_err = sp.stats.linregress(chl_anomalies[:,lat,lon],sst_anomalies[:,lat,lon])
except ZeroDivisionError:
r_value = np.nan
pp = np.nan
#print(lat,lon)
result_cor.append((r_value, float(chl_anomalies.lat[lat]), float(chl_anomalies.lon[lon])))
p_value.append((pp, float(chl_anomalies.lat[lat]), float(chl_anomalies.lon[lon])))

result_cor=pd.DataFrame(result_cor, columns=('correlation', 'lat', 'lon'))
p_value=pd.DataFrame(p_value, columns=('p_value', 'lat', 'lon'))

pv_result_cor= pd.pivot_table(result_cor, index=result_cor.lat, columns=result_cor.lon,values='correlation')
xx,yy = np.meshgrid(anomalies_chl.lon, anomalies_chl.lat)
X1, Y1 = np.meshgrid(pv_result_cor.columns,pv_result_cor.index)

PLOT

fig = plt.figure(figsize=(14,8),constrained_layout=True)
#ax = plt.axes(projection=ccrs.PlateCarree())
ax=plt.axes(projection=ccrs.Mercator(central_longitude=7.1))
img=sst_mean.plot.pcolormesh(cmap='Wistia', shading='auto',add_labels=False,
cbar_kwargs={'shrink': 0.99,'pad':0.01},extend='neither',transform=ccrs.PlateCarree())
ax.set_xlabel('Longitude (°)')
ax.set_ylabel('Latitude (°)')
ax.add_feature(cfeature.LAND, facecolor='grey')
ax.coastlines()
ax.set_yticks([53, 54,55,56], crs=ccrs.PlateCarree()) # -90, -60, -30
ax.set_xticks([3, 6, 9], crs=ccrs.PlateCarree())
ax.set_yticklabels([53, 54,55,56]) # -90, -60, -30
ax.set_xticklabels([3, 6, 9])
fig.text(0.782, 0.89, '(m)')

[Asheef1447]

Thank you so much for taking the time to post this. Would you mind telling me how you read the data in this code? especially this part:

Correlation
anomalies_chl ('time: 276,lat: 288,lon: 648)
anomalies_sst ('time: 276,lat: 288,lon: 648)

Is it a netcdf file, which was read with xarray? I really appreciate your help.

[fllamorim]

YOu are welcome. Yes, the data is netcdf read with xarray dataarray.

[Asheef1447]

Thank you for your suggestion. I’ve implemented the code you provided, making several adjustments to accommodate my dataset. Despite these changes, I’ve encountered an error: TypeError: Dimensions of C (62, 360) should be one smaller than X(361) and Y(173) while using shading='flat'. This error message suggests referring to the help(pcolormesh) documentation.

I would be extremely grateful if you could take the time to review my modified code. Your expertise would be invaluable in helping me identify any potential mistakes I might have made during the adaptation process:

`import matplotlib.pyplot as plt
import numpy as np
import scipy as sp
import pandas as pd
import xarray as xr
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

Define the paths to the NetCDF files

path_ice = "/content/drive/My Drive/Data/extracted_sea_ice.nc"
path_carbon = "/content/drive/My Drive/Data/regrid_carbon.nc"

Loading data into xarray DataArrays

ds1 = xr.open_dataset(path_ice)
ds2 = xr.open_dataset(path_carbon)

Extract the data variables

sic = ds1['SIC']
spco2 = ds2['spco2']

Calculate correlations and p-values

result_cor = []
p_value = []

for lat in range(len(ds1['LAT'])):
for lon in range(len(ds1['LON'])):
try:
corr, p = sp.stats.pearsonr(sic[:, lat, lon], spco2[:, lat, lon])
except ValueError:
corr = np.nan
p = np.nan
result_cor.append((corr, float(ds1['LAT'][lat]), float(ds1['LON'][lon])))
p_value.append((p, float(ds1['LAT'][lat]), float(ds1['LON'][lon])))

Create DataFrames

result_cor_df = pd.DataFrame(result_cor, columns=('correlation', 'lat', 'lon'))
p_value_df = pd.DataFrame(p_value, columns=('p_value', 'lat', 'lon'))

Create a pivot table for correlation

pv_result_cor = pd.pivot_table(result_cor_df, index=result_cor_df['lat'], columns=result_cor_df['lon'], values='correlation')

Create the plot

fig = plt.figure(figsize=(14, 8), constrained_layout=True)
ax = plt.axes(projection=ccrs.PlateCarree())
img = ax.pcolormesh(ds1['LON'], ds1['LAT'], pv_result_cor.values, cmap='Wistia', shading='auto')
ax.set_xlabel('Longitude (°)')
ax.set_ylabel('Latitude (°)')
ax.add_feature(cfeature.LAND, facecolor='grey')
ax.coastlines()
ax.set_yticks([53, 54, 55, 56], crs=ccrs.PlateCarree())
ax.set_xticks([3, 6, 9], crs=ccrs.PlateCarree())
ax.set_yticklabels([53, 54, 55, 56])
ax.set_xticklabels([3, 6, 9])
fig.colorbar(img, ax=ax, shrink=0.99, pad=0.01, extend='neither', label='Correlation')

Show the plot

plt.show()
`

[ahuang11]

Hvplot_obj * hvplot_obj2 should overlay Then hvplot_overlay + hvplot_overlay lays out. Check out holoviews layout (hvplot is method based holoviews)

…

On Tue, Nov 23, 2021, 12:38 AM HARSHBARDHAN KUMAR ***@***.***> wrote: Still I am looking some robust way to include significance plot as dot mark (p<0.05), wherever applicable over spatial trend plot. Suggestion regarding *hvplot.points*, seems not compatible with the grouped sub-plots (e.g., 3row*3col=9 plots) . hvplot can be used with single plot. It also takes enough memory and time to render the single plot. — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#5982 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADU7FFTY2B7EJQKBDSR6FCLUNMZGVANCNFSM5H6KXA2Q> .

[Mohima17]

https://github.com/royalosyin/Python-Practical-Application-on-Climate-Variability-Studies/blob/master/ex34-Correlations%20between%20SOI%20and%20SLP%2C%20Temperature%20and%20Precipitation.ipynb

[fllamorim]

Try to transpose your pv table (pv.T) or do xx,yy=np.meshgrid(d1.LON,d1.LAT) https://stackoverflow.com/questions/24791614/numpy-pcolormesh-typeerror-dimensions-of-c-are-incompatible-with-x-and-or-y

…

On Thu, Sep 7, 2023, 22:44 Md. Ashif Imam Khan ***@***.***> wrote: Thank you for your suggestion. I’ve implemented the code you provided, making several adjustments to accommodate my dataset. Despite these changes, I’ve encountered an error: TypeError: Dimensions of C (62, 360) should be one smaller than X(361) and Y(173) while using shading='flat'. This error message suggests referring to the help(pcolormesh) documentation. I would be extremely grateful if you could take the time to review my modified code. Your expertise would be invaluable in helping me identify any potential mistakes I might have made during the adaptation process: `import matplotlib.pyplot as plt import numpy as np import scipy as sp import pandas as pd import xarray as xr import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter Define the paths to the NetCDF files path_ice = "/content/drive/My Drive/Data/extracted_sea_ice.nc" path_carbon = "/content/drive/My Drive/Data/regrid_carbon.nc" Loading data into xarray DataArrays ds1 = xr.open_dataset(path_ice) ds2 = xr.open_dataset(path_carbon) Extract the data variables sic = ds1['SIC'] spco2 = ds2['spco2'] Calculate correlations and p-values result_cor = [] p_value = [] for lat in range(len(ds1['LAT'])): for lon in range(len(ds1['LON'])): try: corr, p = sp.stats.pearsonr(sic[:, lat, lon], spco2[:, lat, lon]) except ValueError: corr = np.nan p = np.nan result_cor.append((corr, float(ds1['LAT'][lat]), float(ds1['LON'][lon]))) p_value.append((p, float(ds1['LAT'][lat]), float(ds1['LON'][lon]))) Create DataFrames result_cor_df = pd.DataFrame(result_cor, columns=('correlation', 'lat', 'lon')) p_value_df = pd.DataFrame(p_value, columns=('p_value', 'lat', 'lon')) Create a pivot table for correlation pv_result_cor = pd.pivot_table(result_cor_df, index=result_cor_df['lat'], columns=result_cor_df['lon'], values='correlation') Create the plot fig = plt.figure(figsize=(14, 8), constrained_layout=True) ax = plt.axes(projection=ccrs.PlateCarree()) img = ax.pcolormesh(ds1['LON'], ds1['LAT'], pv_result_cor.values, cmap='Wistia', shading='auto') ax.set_xlabel('Longitude (°)') ax.set_ylabel('Latitude (°)') ax.add_feature(cfeature.LAND, facecolor='grey') ax.coastlines() ax.set_yticks([53, 54, 55, 56], crs=ccrs.PlateCarree()) ax.set_xticks([3, 6, 9], crs=ccrs.PlateCarree()) ax.set_yticklabels([53, 54, 55, 56]) ax.set_xticklabels([3, 6, 9]) fig.colorbar(img, ax=ax, shrink=0.99, pad=0.01, extend='neither', label='Correlation') Show the plot plt.show() ` — Reply to this email directly, view it on GitHub <#5982 (reply in thread)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AVADOZVMQYD76E6KDOJOK2TXZIWZBANCNFSM5H6KXA2Q> . You are receiving this because you commented.Message ID: ***@***.***>