Error plotting countries spanning both eastern and western hemispheres

[ritviksahajpal]

Thanks for an excellent package! I found that trying to use salem to plot a map for an individual country fails for a country like Russia which wraps around to exist in both western and eastern hemispheres.

borders = borders.loc[borders['name'].isin(['Russia'])]
arr.salem.subset(shape=borders, margin=2)

[fmaussion]

Thanks for the report! Indeed I haven't looked much how salem behaves across hemispheres. I'll try to have a look soon.

Is your data defined on [0, 360] or [-180, 180]?

[ritviksahajpal]

Thanks for looking into this! The shapefile is define on -180 to 180 (and so is the data). The shapefile is here: https://www.dropbox.com/sh/yp22wdue7ri63z9/AADqm4zPwKNkykcneMwthm8Fa?dl=0

[fmaussion]

OK, so there are a couple of things going on. If everything is in the -180 180 range, roi and subset work fine, as long as one uses no margin:

import xarray as xr
import numpy as np
import salem
from salem import get_demo_file, open_wrf_dataset
import matplotlib.pyplot as plt

borders = salem.read_shapefile(get_demo_file('world_borders.shp'), cached=True)
borders = borders.loc[borders['CNTRY_NAME'].isin(['Russia'])]

# make dummy dataset
lon = np.linspace(1, 359, 180) - 180
lat =  np.linspace(1, 89, 45)
data = lon * np.ones((len(lat), 1))
da = xr.DataArray(data, dims=['lat', 'lon'], coords={'lat': lat, 'lon':lon})

# plot
f = plt.figure(figsize=(12, 4))
roi = da.salem.roi(shape=borders)
roi.salem.quick_map(countries=False)

f = plt.figure(figsize=(12, 4))
s = roi.salem.subset(shape=borders)
s.salem.quick_map(countries=False)

With margin=2 we have an issue, as you noticed. I am not sure what to do in this case, though. The subset (with or without margin) isn't very useful...

The easiest way to make Russia "continuous" is to shift the data prior analysis:

sdata = np.roll(data, 90, axis=1)
slon = np.roll(lon, 90)
slon = np.where(slon < 0, slon + 360, slon)
sda = xr.DataArray(sdata, dims=['lat', 'lon'], coords={'lat': lat, 'lon':slon})

But then, the shapefile also has to be converted too :-(

def trafo_360(x, y, z=None):
    return np.where(x < 0, x+360, x), y
from shapely.ops import transform
sborders = borders.copy()
sborders['geometry'] = borders.geometry.apply(lambda geom: transform(trafo_360, geom))

# finally!
f = plt.figure(figsize=(12, 4))
roi = sda.salem.roi(shape=sborders)
s = roi.salem.subset(shape=sborders, margin=2)
s.salem.quick_map(countries=False)

I guess that salem could take over the last step (the shapefile transformation), and could provide automatic tools to shift datasets (I've needed this several times).

[ritviksahajpal]

thank you so much! this is a very comprehensive answer and works for me.