Wradlib.vis.plot_scan_strategy - wrong altitude

Hello everyone,

I am new to this area, but i am wondering if someone else has as well encountered the error message by plotting the scan_strategy-function?

Plotting the scan strategy with: ax = wradlib.vis.plot_scan_strategy(ranges, elevs, sitecoords, beamwidth=vol.how.beamwH) works fine. But if I want to use:

ax = wradlib.vis.plot_scan_strategy(ranges, elevs, sitecoords, cg=True, terrain=True) → I get: TypeError: download_srtm() takes 2 positional arguments but 3 were given.

Now, ax = wradlib.vis.plot_scan_strategy(ranges, elevs, sitecoords, beamwidth=vol.how.beamwH, cg=True) runs fine BUT the altitude is shifted by -1000 as MSL is now there.

The variables in the scan strategy are given by:

‘’’
ranges = swp_0.range.values
sitecoords = (
swp_0.longitude.values.item(),
swp_0.latitude.values.item(),
swp_0.altitude.values.item(),
)
sweeps = list(vol.keys())[:-3]
elevs =
for dataset in sweeps:
elangle = vol[dataset][“where”].elangle
elevs.append(elangle)
‘’’

I hope that someone can help me out on how to fix this issue.

Best, Julian

This should be fixed in wradlib 2.0.2

But you’ll need a EARTHDATA account correctly setup:

See Plotting Radar Scan Strategy — wradlib.

HTH,
Kai

Hi Kai,

thanks a lot. I updated wradlib to 2.0.2, but I still get some errors. If I use:

ax = wradlib.vis.plot_scan_strategy(ranges,
                                    elevs,
                                    sitecoords,
                                    cg=True,
                                    az=180,                                    
                                    cmap="viridis",
                                    units='km',
                                    terrain=True)

everything works “fine”, but the MSL is still leveled at -1km. WHY? It is the same in the last figure of the wradlib-link you added above, and I don’t understand why it should be like this?

Furthermore, if i change az to something different than az=180, it gives me the following error message:

RuntimeError: No usable source images.

I am happy to hear from you,
Julian

Hi Julian,

not sure, it would be great, if you can make a little MCVE from your example. So at least give the input data you have used. This will help to reproduce issue and find a way around it.

This might be a bug using units='km'. Does it work if you omit that kwarg? Another possible cause might be broken site_coords.

Best,
Kai

That’s what I have so far:

import xarray as xr
import datatree
import wradlib
import os

os.environ["WRADLIB_EARTHDATA_BEARER_TOKEN"] = "my_bearer_token"
os.environ["WRADLIB_DATA"] = "/path/to/my/folder"

fname = "/home/ba9817/phd/finland/radar/hdf5/20220126/202201260930_fiika_PVOL.h5" # can be downloaded here: http://fmi-opendata-radar-volume-hdf5.s3-website-eu-west-1.amazonaws.com/?prefix=2022/01/26/fiika/

vol = datatree.open_datatree(fname)

# Calculating the elevation angles
sweeps = list(vol.keys())[:-3] 
elevs = []
for dataset in sweeps:
    elangle = vol[dataset]["where"].elangle
    elevs.append(elangle)

# Reading data from first sweep
swp_0 = xr.open_dataset(fname, engine="odim")

ranges = swp_0.range.values
sitecoords = (
    swp_0.longitude.values.item(),
    swp_0.latitude.values.item(),
    swp_0.altitude.values.item() 
)

ax = wradlib.vis.plot_scan_strategy(ranges,
                                    elevs,
                                    sitecoords,
                                    maxalt=15.0,
                                    cg=True,
                                    az=180,                                    
                                    cmap="viridis",
                                    units='km',
                                    terrain=True)

@lou Code looks good. The RuntimeError: No usable source images. might it have to do with the coverage of the SRTM data? IIRC it’s only covered to 60deg north.

Do you have a local DEM covering whole of Finland? We might just skip the automated SRTM download and prepare the terrain ourselves from local DEM data. If you need help with that you might have a glimpse at Beam Blockage Calculation using a DEM — wradlib how to interpolate a given DEM
onto the radar polar layout.

The other thing, maxalt would still have to be given in meters (even if you specify ‘km’).

The altitude-problem is probably a bug in the curvelinear grid code, as the issue also shows up in the docs. Good catch! Would you mind opening an issue on the wradlib github issue tracker Issues · wradlib/wradlib · GitHub?

Confirmed, it’s a bug, we get this offset of 1000m because we use two different earth radius in calculation and plotting (6370000 and 6371000). If you feel capable enough you can change the values to 6371000 as a quick fix:

Indeed, I forgot about that the SRTM data covers only the area up to 60deg north. I don’t have a local DEM covering Finland - I’ll have a look if I can find one and implement it.

Good to know about maxalt!

I think I don’t feel capable of fixing this in the code as I don’t know what the right way would be to do so.

@lou The altitude issue is fixed in wradlib 2.0.3, already available in PyPI. Next day in conda-forge.

You might wanna take a look at this topography:
https://portal.opentopography.org/datasets?minX=7.2674560546874565&minY=63.01448255321023&maxX=41.27563476562497&maxY=71.86221497468833

I have as well found: Karttapaikka - Maanmittauslaitos

I downloaded some Geo-Tiff files that I need. But I could not figure out how I should now process these data in order to implement it into the scan_strategy-function? Would you mind to help me out?

And, @kmuehlbauer, is there a way to add a star or a dot on the MSL (e.g. a location that is 60km away from the radar station) with some text?

ax = wradlib.vis.plot_scan_strategy(ranges,
                                    elevs,
                                    sitecoords,
                                    maxalt=15000.0,
                                    cg=True,
                                    az=10,                                    
                                    cmap="viridis",
                                    units='km',
                                    terrain=False)
# get polar parasite axis
paax = ax.parasites[1]
# use this as earth radius
er = 6371000
# the distance in m
rng = 60000
# calculate bin distance
ade = wradlib.georef.bin_distance(rng, 0, sitecoords[2], ke=1.0)
# get angle
thetap = -np.degrees(ade / er) + 90.0
# plot the stuff at the wanted location
paax.plot(thetap, er, "r*", linewidth=3, label="RAD", zorder=3)
paax.text(thetap, er, "Radar", zorder=3)

Takes a while, need to find a bit of time doing that. You might explore GDAL merge gdal_merge.py — GDAL documentation on how to merge the tiles. After that just follow-up with the preprocessing steps Beam Blockage Calculation using a DEM — wradlib how to interpolate the DEM onto the radar polar grid. Then just feed the wanted azimuth to terrain (ax = wrl.vis.plot_scan_strategy(ranges, elevs, site, units="km", terrain=polarvalues[wanted_azi]).

HTH,
Kai

Thanks. That workes perfeclty :slight_smile:

Maybe, you can have a look - when you have a spare minute - as well at:

ax = wradlib.vis.plot_scan_strategy(ranges,
                                    elevs,
                                    sitecoords,
                                    maxrange=120000.0,
                                    cg=True,
                                    az=180,
                                    cmap="viridis",
                                    units='km')

For me, maxrange does not work when I use cg=True!

Thanks @lou, indeed, maxrange is not taken into account currently. Let’s put it on the list.

A post was split to a new topic: Unexpected output with wradlib.vis.plot_scan_strategy