7

I am trying to plot a netcdf file from ECMWF ERA5 in a folium map but I am having troubles with the projections in the file.

ECMWF documentation states: ECMWF data is referenced in the horizontal with respect to the WGS84 ellipse (which defines the major/minor axes) but in the vertical it is referenced to the Geoid (EGM96).

The netcdf file in fact seems to line up well in the horizontal axis but not on the vertical axis.

Is it possible to convert the vertical projection to match the horizontal projection in WGS84?

Here is my python code:

import pandas as pd
from streamlit_folium import folium_static
import matplotlib.pyplot as plt
import folium
import numpy as np
import xarray as xr
import pyproj
from pyproj.crs import CompoundCRS

cmpd_crs = CompoundCRS(name="WGS 84 + EGM96 height", components=["EPSG:4326", "EPSG:5773"])
print(cmpd_crs)

ds = xr.open_dataset(r'/home/giacomo/anaconda3/envs/ERA5plot/2021_daily_max.nc')

ds1 = ds.reduce(np.nansum, 'expver')
print('netcdf opened')
print(ds)

sliced = ds1.sel(time=slice(pd.to_datetime('2021-01-28'), pd.to_datetime('2021-02-02')))
print(sliced)

map_data = sliced['i10fg'].max('time')
print(map_data)

np_image = map_data.values

### projection transformation
map_data.rio.set_crs(cmpd_crs)

inProj = pyproj.Proj(cmpd_crs)
outProj = pyproj.Proj('epsg:4326')

x_1, y_1 = map_data['longitude'], map_data['latitude']
xy = list(zip(x_1, y_1))
x, y = pyproj.transform(inProj, outProj, xy[0], xy[1])

### set coloration
val_min = float(np.amin(np_image))
val_max = float(np.amax(np_image))

cm = plt.cm.get_cmap('viridis')

np_image_norm = (np_image - val_min) / (val_max - val_min)

colored_data = cm(np_image_norm)


### create folium map
map_loc = folium.Map(location=[0, 0], zoom_start=0, crs='EPSG3857')
folium.raster_layers.ImageOverlay(
                name="Mercator projection SW",
                image= colored_data, #np_image,
                bounds=[[float(y_1.min()), float(x_1.min())], [float(y_1.max()), float(x_1.max())]],
                opacity=0.5,
                mercator_project = True,
            ).add_to(map_loc)
folium_static(map_loc)

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.