""" import numpy as np import numpy.ma as ma import xarray as xr import esmlab import cartopy.crs as ccrs from cartopy.mpl.geoaxes import GeoAxes from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import matplotlib.pyplot as plt import prettyplotlib as ppl from prettyplotlib import brewer2mpl from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection from mpl_toolkits.basemap import shiftgrid from mpl_toolkits.basemap import Basemap import matplotlib.colors as mcolors from mpl_toolkits.axes_grid1 import make_axes_locatable from cartopy.util import add_cyclic_point from scipy.stats.stats import pearsonr import scipy.stats as stats from esmtools.stats import* import matplotlib.patches as mpatches """ ds3 = xr.open_dataset('/tigress/sirishak/postpdata/urban_proj/hist/urban/land_mask.land_mask.nc') land_mask = ds3['land_mask'].sel(lon=slice(360-130, 360-60), lat=slice(20, 55)) #land_mask = ds3['land_mask'].sel(lon=slice(360-125, 360-117), lat=slice(43, 54)) #years = slice(1990, None) def is_amj(month): return (month >= 6) & (month <= 8) #return (month >= 3) & (month <= 5) #return (month < 3) | (month == 12) # DJFM #ds1 = xr.open_dataset('/tigress/sirishak/postpdata/urban_proj/hist/urban/PNW_data/noboundary_flux/maps/t.urban_PNW_combined_no_boundary_flux_NA_land.nc') ds1 = xr.open_dataset('/tigress/sirishak/postpdata/urban_proj/hist/urban/PNW_data/correct_data/wastecool_global_urbanisation_00_11_14_all_ens/t_ref.urban_PNW_1990_2014_global_urbanisation_all_ens.nc') ds1 = ds1['t_ref'].sel(lon=slice(360-130, 360-60), lat=slice(20, 55))#*86400 ds1=ds1.sel(time=slice('1990-01-01', '2014-12-31')) #ds1 = ds1.sel(lon=slice(360-125, 360-117), lat=slice(43, 54))#.sel(lon=slice(360-130, 360-60), lat=slice(20, 55))#.sel(time=slice('1990-01-01', '2020-01-01')) #ds1= ds1.groupby('time.year').max('time')#.sel(time=slice('1980-01-01', '2010-01-01'))#.mean(dim='time') ds1 = ds1.sel(time=is_amj(ds1['time.month'])).groupby('time.year').mean('time') #clint(ds1) ds1.coords['land_mask'] = (('lat', 'lon'), land_mask.data) print(ds1) ds1 = ds1.where(ds1.land_mask > 0.3) ds1 = ds1.stack(s=['en', 'year']) #ds3 = ds1.mean('lon') print(ds1) ds2 = xr.open_dataset('/tigress/sirishak/postpdata/urban_proj/hist/urban/PNW_data/correct_data/nourban_00_21_24_all_ens/t_ref.nourban_PNW_1990_2014_all_ens.nc') ds2 = ds2['t_ref'].sel(lon=slice(360-130, 360-60), lat=slice(20, 55)) ds2=ds2.sel(time=slice('1990-01-01', '2014-12-31')) ds2 = ds2.isel(time=is_amj(ds2['time.month'])).groupby('time.year').mean('time')#.mean('time') ds2.coords['land_mask'] = (('lat', 'lon'), land_mask.data) ds2 = ds2.stack(s=['en', 'year']) print(ds2) ds2 = ds2.where(ds2.land_mask > 0.3) print(ds2) #ds4 = ds2.mean('lon') #levels=[-4,-3,-2,-1,0,1,2,3,4] #levels=[0,20,40,60,80,100] #levels=[-6,-4,-2,-1,0,1,2,4,6] #=[-0.6,-0.4,-0.2,-0.1,0,0.1,0.2,0.4,0.6] #levels=[-10,-8,-6,-4,-2,0,2,4,6,8,10] #levels=[-1.5,-1,-0.5,-0.25,0,0.25,0.5,1,1.5] levels=[-0.4,-0.3,-0.2,-0.1,0,0.1,0.2,0.3,0.4] da = (ds1 - ds2)#*10E-3#*2.5E6 sys.exit() #da = (ds1 - ds3) - (ds2 - ds4) #/((ds1+ds2)/2))*100) #*2.5E6 print(da.shape) #da = da.mean('en') red_blue =brewer2mpl.get_map('RdBu', 'diverging', 11,reverse=True).mpl_colormap red_y_green = brewer2mpl.get_map('BrBG', 'diverging', 11).mpl_colormap PuOr=brewer2mpl.get_map('PuOr', 'diverging', 11).mpl_colormap YlOrBr=brewer2mpl.get_map('YlOrBr', 'sequential', 9).mpl_colormap cmap =red_blue norm = mcolors.BoundaryNorm(levels, cmap.N, clip=True) #extent = mtransforms.Bbox.union([col.get_datalim(self.transData) # for col in result.collections]) from scipy.stats import sem, t from scipy import mean confidence = 0.90 n = da.s.size m = da.mean('s') std_err = sem(da, axis=2) h = std_err * t.ppf((1 + confidence) / 2, n - 1) print(m) da_h = xr.DataArray(h, dims=('lat', 'lon'), coords=(da.lat, da.lon)) proj = ccrs.PlateCarree() #proj = ccrs.Orthographic(40, 90) mask = np.abs(m)>np.abs(da_h) #fig, ax = plt.subplots(nrows=1,subplot_kw=dict(projection=ccrs.NorthPolarStereo())) fig, ax = plt.subplots(nrows=1,subplot_kw=dict(projection=proj)) h = m.plot.contourf(ax=ax, transform=proj, cmap=cmap,levels=levels,extend='both',add_colorbar=False) m.where(mask).plot.contourf(ax=ax, transform=proj, colors='none', hatches=["...."], add_colorbar=False) #plt.colorbar() #divider = make_axes_locatable(ax) #ax_cb = divider.new_horizontal(size="3%", pad=0.3, axes_class=plt.Axes) #ax.set_xticks(np.linspace(-180, 180, 13), crs=proj) #ax.set_yticks(np.linspace(-90, 90, 11), crs=proj) ax.set_xticks(np.linspace(-130, -60, 8), crs=proj) ax.set_yticks(np.linspace(20, 55, 8), crs=proj) lon_formatter = LongitudeFormatter() lat_formatter = LatitudeFormatter() ax.xaxis.set_major_formatter(lon_formatter) ax.yaxis.set_major_formatter(lat_formatter) ax.gridlines(color='gray', alpha=0.6, linestyle='--') ax.coastlines() #ax.set_extent([-180, 180, 20, 90], crs=ccrs.PlateCarree()) #ax.set_extent([-180, 180, -90, 90]) #fig.add_axes(ax_cb) #plt.colorbar(h, cax=ax_cb) fig.colorbar(h, orientation='horizontal', pad=0.08) patches = [] zone_A = np.array([[-125,54],[-117,54],[-117,43],[-125,43]]) patches.append(Polygon(zone_A)) ax.add_collection(PatchCollection(patches, facecolor='None', edgecolor='k', linewidths=1.0)) plt.tight_layout() plt.show()