# Extracting Regional Data from ERA5

This is coming from a newcomer, so please be kind.

As I delve into the world of meteo data and reanalysis for the first time, I found that acquiring weather data is a lengthy process, so I'm shouting out for some help to get me started with it.

1. What is the recommended workflow to acquire ERA5 data for a region of interest (ROI) for a few surface variables (precip, soilmoisutre, tempmin, solar rad)?

• Is there an API that connects me (R, Python, Curl) to the ERA5 database?
• Do I need to download whole (global) files even if my ROI is always within the same country? If so, are we talking about 100GB, 1,000GB or 10,000GB of data for 6 variables, whole historical?
2. Seems like there are plenty of tools to extract ROI data from grib files. What is the recommended tool to learn in 2019?

• Is CDO the most appropriate tool for the job?

• Any Python/R packages that can extract and plot these data in an efficient way?

Thanks and let me know how to complement the info above

• Do you need to automate the process? Or you are happy with an online interface to manually select, visualize and download the data? – Camilo Rada Mar 26 at 17:16
• if it is just few layers for a specific area, you can simply follow the procedure described here (confluence.ecmwf.int/display/CKB/How+to+download+ERA5) and use the built in API. – Nemesi Mar 26 at 17:19
• @CamiloRada, data acquisition has to be automated. I've found the manual interface <<< confluence.ecmwf.int/display/CKB/How+to+download+ERA5>> and studying it now. Question remains on what the best tool to aggregate this data for ROIs is. Would you suggest ClimateDataOperators (CDOs) for someone on Windows ? – Dan Mar 27 at 17:44
• @Dan, Sorry, I didn't read all the tutorial and did not realize that the facility was still under development. If you are familiar with Python, you could adapt this script (confluence.ecmwf.int/pages/viewpage.action?pageId=73016800) to extract the data for your ROI. I did not not work yet with ERA5 files, but I assume that the structure of the file should be similar to the one they used for ERA-Interim. – Nemesi Mar 28 at 9:36
• To answer your question, in my experience, CDO and NCO are the most efficient tools to resize and get info from NetCDF, HDF5, or GRiB files, but sometimes you want to use scripts (I usually use R, but Python is probably even more efficient) to extract and process climate data and prepare the outputs in a format that you can use for your purposes. For instance, I often used R scripts to extract climate projections from .nc files and prepare input files for a land surface model. So, there is no a standard way to proceed, it depends on what you need to do. – Nemesi Mar 28 at 9:40

To me I am not sure why the subset option should not work with a python script. According to me it does work. Here is a sample script.

import cdsapi

c = cdsapi.Client()

c.retrieve(
'reanalysis-era5-single-levels',
{
'variable':'total_precipitation',
'product_type':'reanalysis',
'year':'2010',
'month':'04',
'day':'07',
'area':'60/0/0/100',
'time':[
'00:00','01:00','02:00',
'03:00','04:00','05:00',
'06:00','07:00','08:00',
'09:00','10:00','11:00',
'12:00','13:00','14:00',
'15:00','16:00','17:00',
'18:00','19:00','20:00',
'21:00','22:00','23:00'
],
'format':'netcdf'
},
'precip.nc')


UPDATE In response to some of the comments raised by @Nemesi - with this CDS Python API There is no need to use CDO or NCO to subset the data. The API itself takes care of it for you using the keyword area. Once you specify the bounding box the data that goes get downloaded is the bounding box boundaries you specify. There is also NO programmatic subsetting of global data. Right now the CDS web interface does not allow you to do subsetting(the problem is being fixed as I write according to Copernicus support). But that is not really a problem. Since OP's requirement is one of automating the process the CDS Python Web API fills the purpose. In this context subset is the same as ROI or region of interest.

Here is what I had to do get this to work. I had to install the cdsapi for python. One can do this using pip (I usually do it - python3.6 setup.py install since I do not have conda). Then create a .cdsapirc(on Linux and other UNIX flavors it should be under the main home directory. Create this file wherever you have the \$HOME variable defined) which should be like this -

url: https://cds.climate.copernicus.eu/api/v2
key: {UID}:{API key}
verify:0


The values of UID and API key you should get when you register at the new CDS web interface site.

And then just run it.

Here the key parameter is the area

and the first two values are starting latitude and longitude followed by the ending latitude and longitude. When I viewed the downloaded netCDF file using ncdump it did have the right bounding box values. Now you could change it to GRIB format if you so wish but I think the basic format of the script does not change.

Here is a plot of the subset of the precipitation for a specific instant in a 24 hour period.

Here is the code that plots it (I use matplotlib and cartopy)

from netCDF4 import Dataset,num2date
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib as mpl
mpl.rcParams['mathtext.default'] = 'regular'
import matplotlib.pyplot as plt
import numpy as np

file = "precip.nc"
nc_pvFile = Dataset(file,'r')
lats = nc_pvFile.variables['latitude'][:]  # extract/copy the data
lons = nc_pvFile.variables['longitude'][:]
lats = lats[:].squeeze()
lons = lons[:].squeeze()

preciPlot = nc_pvFile.variables['tp'][:]

pp = preciPlot[0,:,:]

ax1 = plt.axes(projection=ccrs.PlateCarree(central_longitude=180))
clevs = np.arange(min(pp.flatten()),max(pp.flatten())*1000,1)

shear_fill = ax1.contourf(lons,lats,pp*1000,clevs,
transform=ccrs.PlateCarree(),cmap=plt.get_cmap('hsv'),
linewidth=(10,),levels=100,extend='both')

ax1.coastlines()
ax1.gridlines()
ax1.set_xticks([0, 10,20,30,40,50,60,70,80,90,100], crs=ccrs.PlateCarree())
ax1.set_yticks([0, 10,20,30,40,50,60], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True,
number_format='.0f')
lat_formatter = LatitudeFormatter()
ax1.xaxis.set_major_formatter(lon_formatter)
ax1.yaxis.set_major_formatter(lat_formatter)
cbar = plt.colorbar(shear_fill, orientation='horizontal')
plt.title('Total precipitation', fontsize=16)
plt.savefig('precip_era.png')


• Hi @gansub, no one said that Python (or R or Matlab) are not able to subset a GRIB or a NetCDF (I actually use R all the times to do that), I just said that NCO and CDO are more efficient. – Nemesi Apr 2 at 7:28
• @Nemesi - No. Please read the question again. The CDS web interface for download of ERA interim data does NOT do a subset currently(the problem is being fixed according to ECMWF support). It only allows full global data downloads. My answer says NO that is incorrect. Using the area keyword you can subset using the automated process. If this premise is incorrect I will delete my answer. I will let OP clarify that in a comment or an edit to his question. – gansub Apr 2 at 8:27
• @Nemesi It is not about Python or R doing the subsetting. The CDS API allows you to subset without using programmatic functionality. – gansub Apr 2 at 8:32
• don't get me wrong, yours is a good answer and I up-voted it. I am actually not familiar with the CDS API, and I haven't yet used the ERA5 data. Glad to hear that now works. Thanks for address this issue. – Nemesi Apr 2 at 8:55
• @Dan There is NO local extraction of the data. It is completely remote. Please install CDS python API and come to the conclusion yourself. – gansub Apr 2 at 13:58

Any Python/R packages that can extract and plot these data in an efficient way?

I can highly recommend using the Python package xarray for data analysis and Cartopy for plotting.