plotting ppi on a basemap (not! using the "basemap.contourfill" function)

838 views
Skip to first unread message

Dzengiz Tafa

unread,
Mar 27, 2014, 12:57:38 AM3/27/14
to wradli...@googlegroups.com
Hi everybody.

I am new to the group & wradlib so excuse me for a possible stupid question.
I find myself struggling to plot a ppi on a basemap. eg a ppi on a background of the Benelux like the one here (hosted & created by the Bonn Meteorological institute).

I found a way to plot the data from the Jabbeke site in Belgium but only using the "basemap.contourfill" function.
Although it is a good way to show precipitation as a whole I loose a lot of detail. To that extent that I loose the gate to gate signatures on which the presentation I am to give on the 11th of may is based (tornado event of november 3d, 2013). An example of the product using the basemap.contourfill function is here:

Can anyone point me in the right direction as to how I can plot the radar data which gives a presentation like the upper picture (Bonn Meteorological institute)? I tried following the examples but as good as they are, they lead me nowhere.

If somebody could find the time to set me on the right track or perhaps share some code, I would be very grateful as it will surely help me present the best material on the conference I need to attend for my talk about the severeweather outbreak.

Grtz,


Dzengiz

Kai Muehlbauer

unread,
Mar 27, 2014, 3:25:32 AM3/27/14
to wradli...@googlegroups.com
Hi Dzengiz,

you are very welcome.

Your email to radar...@uni-bonn.de reached me too, but I'll answer
here, to cover this for all who are interested in this topic.

As you may have noticed by looking at the MIUB website, we have lots of
radar images for several scans every 5 minutes. Since Basemap is "very"
slow in high resolution mode, I had to take some detours to get things
faster. And of course, I had problems with the correct projection, never
knowing if I got it right in the end.

Here is how I have done it, but I do this on a server without graphical
output using "mpl.use('Agg')" which means plotting only to buffer:

1. create overlay separately:
1a. Plot your radar data in the projection you want, preserve the
axis-object (you could just create the axis and set the xlim and ylim
coordinates accordingly)
1b. Create the basemap according to your projection and declare width
and height correctly
m = Basemap(resolution='h', projection='yourproj', ax=yourax,
width=yourwidth, height=yourheight, lat_0=lat_0, lon_0=lon_0)
1c. plot all overlay you need with the basemap functions (rivers,
countries etc.)
1d. grab the whole figure PIL image as explained here:
http://www.icare.univ-lille1.fr/wiki/index.php/How_to_convert_a_matplotlib_figure_to_a_numpy_array_or_a_PIL_image
1e. crop your image to your axes bounding box, you get only the axis
image then
1f. set the background color to alpha with this nice snippet from Nadia
Alramli:
http://stackoverflow.com/questions/1616767/pil-best-way-to-replace-color
1g. save this image

2. create underlay (orography) separately
2a. here you can create your basemap all the same but you can use
"crude" resolution because we draw the orography from the srtm data
2b. get your srtm data, eg. from here:
http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp
2c. read the srtm data using this snippet here:
https://stevendkay.wordpress.com/tag/srtm/
2d. since you read kind of tiles from the srtm files you have to cut
your viewing part out of the srtm data
2e. create linear spaced lat,lon arrays and plot the data with this
coordinates on the basemap with pcolormesh, I use cm.gist_earth as colormap
2f. repeat steps 1d,e and g

3. plot data with underlay and overlay
3a. use PIL to open your underlay and overlay images
3b. plot your data, preserve axis object
3c. plot underlay with imshow on your axis (zorder=0)
3d. plot overlay with imshow on your axis (zorder=2)
3e. repeat steps 1d and g

I hope you get at least an impression on where to look and how to proceed.

For all: If there is enough interest in that, I could think about
cleaning up my code and put it somewhere in the examples section.

Cheers,
Kai

Am 27.03.2014 05:57, schrieb Dzengiz Tafa:
> Hi everybody.
>
> I am new to the group & wradlib so excuse me for a possible stupid question.
> I find myself struggling to plot a ppi on a basemap. eg a ppi on a
> background of the Benelux like the one here (hosted & created by the Bonn
> Meteorological institute).
>
> <https://lh5.googleusercontent.com/-qiY1MnI2qmo/UzOsjkb7GbI/AAAAAAAAAi4/esOnXNNyxwM/s1600/bonn+radarimage+ppi.png>
>
> I found a way to plot the data from the Jabbeke site in Belgium but only
> using the "basemap.contourfill" function.
> Although it is a good way to show precipitation as a whole I loose a lot of
> detail. To that extent that I loose the gate to gate signatures on which
> the presentation I am to give on the 11th of may is based (tornado event of
> november 3d, 2013). An example of the product using the basemap.contourfill
> function is here:
>
> <https://lh4.googleusercontent.com/-PmUCaBMND6U/UzOuDSBjFjI/AAAAAAAAAjE/udxvl7lRJ-I/s1600/nexrad+kleurenschaal+1dbz+interval.png>
>
> Can anyone point me in the right direction as to how I can plot the radar
> data which gives a presentation like the upper picture (Bonn Meteorological
> institute)? I tried following the examples but as good as they are, they
> lead me nowhere.
>
> If somebody could find the time to set me on the right track or perhaps
> share some code, I would be very grateful as it will surely help me present
> the best material on the conference I need to attend for my talk about the
> severeweather outbreak.
>
> Grtz,
>
>
> Dzengiz
>

--
Kai Muehlbauer
Meteorological Institute University of Bonn
Auf dem Huegel 20 | +49 228 739083
D-53121 Bonn | kai.mue...@uni-bonn.de

Thomas Pfaff

unread,
Mar 27, 2014, 5:21:11 AM3/27/14
to wradli...@googlegroups.com
Hi Kai,

first of all: I don't like basemap, because it's so large and all of it just for the hires coastlines - so I'm a bit biased.

Remembering our discussion yesterday, I also get a bad gut feeling about how basemap might handle projections...

The main benefits of basemap are:
- you don't have to care about the projection
- you have coast- and country borders
- it may make plotting shapefiles easier
The drawbacks:
- it is almost prohibitively large and slow.
- all input must be in lat/lon (and we don't know exactly, if these are WGS84 lat/lons or the dreaded pyproj default)
- the basemap axes object behaves differently from normal axes object, and I believe it does not implement all plotting functions. When I played around with it the last time, manipulating these objects was difficult if not impossible.

If all your data is in the same projection, you don't need basemap and can just plot away with the matplotlib functions.
If you have shapefiles (or another source of your borders, rivers, cities, etc), you can read those using pyshp and then plot that into matplotlib.
After that you have standard matplotlib axes objects, which are nicer to manipulate than the basemap axes, I believe.

There are a few pitfalls especially when using polygon shapes with pyhsp and matplotlib, but these, I think, I've overcome.
Would it be of interest for you and all wradlib users, if we added a lightweight way to produce GIS-like plots? (of course it would be, but would you also be willing to invest the effort :-)).
Or am I totally wrong about all this?

Thanks for the hint about reading SRTM files.

Cheers


Thomas
--
You received this message because you are subscribed to the Google Groups "wradlib-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to wradlib-user...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.



Kai Muehlbauer

unread,
Mar 27, 2014, 5:48:45 AM3/27/14
to wradli...@googlegroups.com, wradl...@googlegroups.com
Please have further discussion on developer list.

Hi Thomas,

your remarks all make sense to me.

I must admit, that I only use basemap for plotting the rivers and
borders. So I could abandon basemap completely when I read that data as
I have done with the srtm data.

In my opinion it would be good if we have at least something which shows
howto produce gis-like plots. I think a joint effort would do the trick.

Since wradlib will rely on gdal/ogr/osr in the future (if Edouards PR is
merged) we also can use all geometry (shapes etc.) from it. For reading,
transforming and plotting shapes you only need "ogr, osr" and
matplotlib.patches "polygon". So I could abandon all extra modules like
shapely, pyshp etc. pp.

What do you think?

Cheers,
Kai

Dzengiz Tafa

unread,
Mar 27, 2014, 1:13:52 PM3/27/14
to wradli...@googlegroups.com

Hi mr Muehlbauer

Thank you so much for your fast reply. I appreciate it a lot.

I fiddled around some more and I came up with this result. However... I guess now all I have to do is make the white pixels transparent?
It sure looks like I'm heading the right way...

Is it possible to confirm that I am doing the right steps?

#basemap
m = Basemap(width=1200000,height=900000,projection='lcc', resolution='i',lat_0=51.182972,lon_0=3.094839)
m.drawcountries(linewidth=1, linestyle='solid', color='black', antialiased=1, ax=None, zorder=1)
m.drawcoastlines(linewidth=1, linestyle='solid', color='black', antialiased=1, ax=None, zorder=2)

ppi = wrl.vis.cg_plot(ind='PPI',fig=fig0,alpha=0.0)
circle = ppi.plot(masked_array, data_range=[0.,180.], theta_range=[0,360], x_range=[-180,180.], x_res = 40., y_range=[-180,180], y_res = 40., radial_range = [0.,180.], a_res = 10, aspect=1.0, faxis=None, alpha=1)

If I should place the code which makes the white pixels transparent as you mentioned (http://stackoverflow.com/questions/1616767/pil-best-way-to-replace-color) I'm guessing it will work?

I am planning to contintue through the course of the evening.

Again:Thank you for you replies guys.

Grtz,


Dzengiz.

Kai Mühlbauer

unread,
Mar 27, 2014, 1:27:11 PM3/27/14
to wradli...@googlegroups.com
Hi Dzengiz,

you can call me Kai, it would make me feel way younger then ;-)

Yes for the first impression it seems that you are on the right track.
The biggest problem is to find the right projection for basemap and your
data.

But anyway, as Thomas pointed out today earlier, basemap is a mess. If
there is anything we can do, just ask around once more. If I find time I
will prepare several code snippets which will handle all those points
and plot in the gui. But this will take a little time.

Cheers,
Kai


Am 27.03.2014 18:13, schrieb Dzengiz Tafa:
> Hi mr Muehlbauer
>
> Thank you so much for your fast reply. I appreciate it a lot.
>
> I fiddled around some more and I came up with this result. However... I
> guess now all I have to do is make the white pixels transparent?
> It sure looks like I'm heading the right way...
>
> <https://lh4.googleusercontent.com/-mzaoKa17Y74/UzRaf0OQzQI/AAAAAAAAAjc/MAWHZX3towg/s1600/closer.png>
>
> Is it possible to confirm that I am doing the right steps?
>
> #basemap
> m = Basemap(width=1200000,height=900000,projection='lcc',
> resolution='i',lat_0=51.182972,lon_0=3.094839)
> m.drawcountries(linewidth=1, linestyle='solid', color='black',
> antialiased=1, ax=None, zorder=1)
> m.drawcoastlines(linewidth=1, linestyle='solid', color='black',
> antialiased=1, ax=None, zorder=2)
>
> ppi = wrl.vis.cg_plot(ind='PPI',fig=fig0,alpha=0.0)
> circle = ppi.plot(masked_array, data_range=[0.,180.],
> theta_range=[0,360], x_range=[-180,180.], x_res = 40.,
> y_range=[-180,180], y_res = 40., radial_range = [0.,180.], a_res = 10,
> aspect=1.0, faxis=None, alpha=1)
>
> If I should place the code which makes the white pixels transparent as
> you mentioned
> (http://stackoverflow.com/questions/1616767/pil-best-way-to-replace-color <http://stackoverflow.com/questions/1616767/pil-best-way-to-replace-color1g>)
> I'm guessing it will work?
>
> I am planning to contintue through the course of the evening.
>
> Again:Thank you for you replies guys.
>
> Grtz,
>
>
> Dzengiz.
>
> On Thursday, March 27, 2014 5:57:38 AM UTC+1, Dzengiz Tafa wrote:
>
> Hi everybody.
>
> I am new to the group & wradlib so excuse me for a possible stupid
> question.
> I find myself struggling to plot a ppi on a basemap. eg a ppi on a
> background of the Benelux like the one here (hosted & created by the
> Bonn Meteorological institute).
>
> <https://lh5.googleusercontent.com/-qiY1MnI2qmo/UzOsjkb7GbI/AAAAAAAAAi4/esOnXNNyxwM/s1600/bonn+radarimage+ppi.png>
>
> I found a way to plot the data from the Jabbeke site in Belgium but
> only using the "basemap.contourfill" function.
> Although it is a good way to show precipitation as a whole I loose a
> lot of detail. To that extent that I loose the gate to gate
> signatures on which the presentation I am to give on the 11th of may
> is based (tornado event of november 3d, 2013). An example of the
> product using the basemap.contourfill function is here:
>
> <https://lh4.googleusercontent.com/-PmUCaBMND6U/UzOuDSBjFjI/AAAAAAAAAjE/udxvl7lRJ-I/s1600/nexrad+kleurenschaal+1dbz+interval.png>
>
> Can anyone point me in the right direction as to how I can plot the
> radar data which gives a presentation like the upper picture (Bonn
> Meteorological institute)? I tried following the examples but as
> good as they are, they lead me nowhere.
>
> If somebody could find the time to set me on the right track or
> perhaps share some code, I would be very grateful as it will surely
> help me present the best material on the conference I need to attend
> for my talk about the severeweather outbreak.
>
> Grtz,
>
>
> Dzengiz
>
> --
> You received this message because you are subscribed to the Google
> Groups "wradlib-users" group.
> To unsubscribe from this group and stop receiving emails from it, send
> an email to wradlib-user...@googlegroups.com
> <mailto:wradlib-user...@googlegroups.com>.
Message has been deleted
Message has been deleted

Kai Muehlbauer

unread,
Mar 28, 2014, 4:24:32 AM3/28/14
to wradli...@googlegroups.com
Hi Dzengiz,

just now I saw that you are trying to use a deprecated function.

The colorbar, colormap and title you get with:

# be aware that the first return is axis object and second is mappable

yourcmap = mpl.cm.jet # use cm according to your needs

ax, ppi = wradlib.vis.plot_ppi(masked_array, r, az, autoext=None,
site=(x,y), elev=0.4, vmin=-0.0, vmax=80.0, cmap=yourcmap)

t = plt.title('Your special Title')
t.set_y(1.05) # adjust spacing as necessary
cbar = plt.gcf().colorbar(ppi, pad=0.075) # adjust padding as necessary
cbar.set_label('reflectivity [dBZ]')
plt.tight_layout() # test tight-layout, it may help beautifying the plot

There are many more possibilities to add colorbars, check also the
CG_Plot tutorial.

Cheers,
Kai




Am 28.03.2014 05:00, schrieb Dzengiz Tafa:
> Finally some results :-)
>
> I made it work on a simple basemap of europe, along with the blue marble
> display. The only thing i can not get to work is my colorbar.
> It seems I have some more fiddling around to do. Here you can see what i
> created. Reflectivity (elevation angle 0.4).
>
> <https://lh5.googleusercontent.com/-FvKPkQOEEY8/UzTlRILkjpI/AAAAAAAAAjs/3sa37twtbnA/s1600/resultaat.png>
>
> As you can see the plot now includes the gates instead of a contourfill of
> the precipitation. Im quite impressed with the result.
>
> However... Is it possible to provide some guidelines as to how I can insert
> a colorbar & title?
>
>
> When erasing the comment char before the colorbar rule (see attachment) I
> get the error message that there is no "autoscale=true" anywhere. I have
> tried to find what it meant, but It seems i can not get it to work.
>
> I would also like to add a title, but when adding a title (or trying to add
> one) i get the error message that "the tuple has no attribute 'title'". I
> tried using the functions in the "vis.py" to plot my scales, title,
> subtitle & date, time but I guess it's because of the way i built up the
> vis.plot_ppi method. Am i correct? If so, is there anyway around it?
>
> Also: when i try to change the colormap to gist_ncar, nothing seems to
> happen. It still keeps its colormap (i think its "jet" but I am not sure of
> this... It Iooks like it though). Is it possible to point me to the
> method/py where the colormap is fixed? I can then easily replace the map
> with the one I want, or just copy the function with another colormap in it.
> No biggy. I like to keep some diverse colormaps, so i can then call the
> right plot-function whenever i feel another colormap is a good idea.
>
> The thing is that I have velocity plots and i can not get the "seismic" map
> to work (you know the red-blue one) on my velocity plots since i dont know
> where i can choose the colormap: as you can see here it's only this blue
> color, no matter what colormap I try to set.
>
> <https://lh6.googleusercontent.com/-vi6OlCDJE0A/UzTxzmmvDJI/AAAAAAAAAj8/1IpTAqCbdS8/s1600/velocity.png>
>
> So the main search now is the colormaps, colorbars & titles. It could very
> well be that i am placing the colormap assignment, colorbar & title wrong,
> but then again: i get no errormessage so... they might be overwritten
> somewhere. If you don't mind i placed an attachment so you can see what i
> did.
>
> Again, thank you very much for your kind help, It has brought me further
> that i'd dare have hoped. It was a rough ride, but now the foundations were
> set, it turned out to be quite smooth sailing.
>
> Grtz
>
> Dzengiz
>
>
> On Thursday, March 27, 2014 5:57:38 AM UTC+1, Dzengiz Tafa wrote:
>>
>> Hi everybody.
>>
>> I am new to the group & wradlib so excuse me for a possible stupid
>> question.
>> I find myself struggling to plot a ppi on a basemap. eg a ppi on a
>> background of the Benelux like the one here (hosted & created by the Bonn
>> Meteorological institute).
>>
>>
>> <https://lh5.googleusercontent.com/-qiY1MnI2qmo/UzOsjkb7GbI/AAAAAAAAAi4/esOnXNNyxwM/s1600/bonn+radarimage+ppi.png>
>>
>> I found a way to plot the data from the Jabbeke site in Belgium but only
>> using the "basemap.contourfill" function.
>> Although it is a good way to show precipitation as a whole I loose a lot
>> of detail. To that extent that I loose the gate to gate signatures on which
>> the presentation I am to give on the 11th of may is based (tornado event of
>> november 3d, 2013). An example of the product using the basemap.contourfill
>> function is here:
>>
>>
>> <https://lh4.googleusercontent.com/-PmUCaBMND6U/UzOuDSBjFjI/AAAAAAAAAjE/udxvl7lRJ-I/s1600/nexrad+kleurenschaal+1dbz+interval.png>
>>
>> Can anyone point me in the right direction as to how I can plot the radar
>> data which gives a presentation like the upper picture (Bonn Meteorological
>> institute)? I tried following the examples but as good as they are, they
>> lead me nowhere.
>>
>> If somebody could find the time to set me on the right track or perhaps
>> share some code, I would be very grateful as it will surely help me present
>> the best material on the conference I need to attend for my talk about the
>> severeweather outbreak.
>>
>> Grtz,
>>
>>
>> Dzengiz
>>
>

Reply all
Reply to author
Forward
0 new messages