I am trying to draw climatology map of GPCC precipitation.
But I found out there is a blank (data cut) in Prime Meridian Line.
How can I fix the problem? I can use CDO, NCO, Python.
I also share the code and data.
(data) https://drive.google.com/drive/folders/1rEHKz5GQlvC3m_Cfzx48cTrPZPU0qAar?usp=sharing
(GPCC metadata) https://opendata.dwd.de/climate_environment/GPCC/html/fulldata-monthly_v2022_doi_download.html
type here
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import calendar
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
import cftime
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
%matplotlib inline
mpl.rcParams['figure.figsize'] = [8., 6.]
filename = 'D:/ERA5/precip.nc'
ds = xr.open_dataset(filename)
ds
da = ds['precip']
da
def is_jjas(month):
return (month >= 6) & (month <= 9)
dd = da.sel(time=is_jjas(da['time.month']))
def is_1982(year):
return (year> 1981)
dn = dd.sel(time=is_1982(dd['time.year']))
dn
JJAS= dn.groupby('time.year').mean('time')
JJAS2 = JJAS.mean(dim='year', keep_attrs=True)
JJAS2
fig, ax = plt.subplots(1, 1, figsize = (16, 8), subplot_kw={'projection': ccrs.PlateCarree()})
cs = plt.contourf(JJAS2.lon, JJAS2.lat, JJAS2, levels=[0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500],
vmin=0, vmax=500, cmap='Blues', extend='both')
# Set the figure title, add lat/lon grid and coastlines
ax.set_title('', fontsize=16)
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
ax.coastlines(color='black')
ax.set_extent([-20, 30, 0, 30], crs=ccrs.PlateCarree())
cbar = plt.colorbar(cs,fraction=0.05, pad=0.04, extend='both', orientation='horizontal')
I tried to search Google and find CDO methods. If I interpolate the data, its resolution might be changed. I want to use this GPCC data with same resolution but no any blanks.
A useful function is add_cyclic
which adds a repeated point for 0 -> 360 longitude. Try something like this:
import cartopy.util as cutil
lon2d, lat2d = np.meshgrid(JJAS2.lon,JJAS2.lat)
cdata, clon2d, clat2d = cutil.add_cyclic(JJAS2, lon2d, lat2d)
And then use these for the plot...