Plotting radar reflectivities in lat lon map

[Krishna chandramouli]

Hi,

I have been trying to plot the radar reflectivities read from Gematronicks radar .vol data.

Thanks to Kai for his help succeeding it.

I have referred some manual and i manually calculated the lat/lon/height corresponding to each reflectivity data from azimuth, elevation and range bins.

Here is the code i used .

I need to plot the data in the same way as the function vis.plot_max_plan_and_vert but with lat and lon ordinates

import wradlib as wrl
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
from osgeo import osr

filename = '/home/krishna/Downloads/IITM_DATA/rawdata/2015120123505200dBZ.vol'
rbdict = wrl.io.read_Rainbow(filename)
lon = rbdict['volume']['sensorinfo']['lon']
lat = rbdict['volume']['sensorinfo']['lat']
sitecoords = (float(lat), float(lon), 53.)

elevs = []
numele = int(rbdict['volume']['scan']['pargroup']['numele'])
slices = rbdict['volume']['scan']['slice']

proj = osr.SpatialReference()
proj.ImportFromEPSG(24383)


for nofslice in range(numele):

    # there are slices with 361 rays
    # make sure to only get first 360 rays
    azi = slices[nofslice]['slicedata']['rayinfo']['data'][0:360]

    azidepth = float(
        slices[nofslice]['slicedata']['rayinfo']['@depth'])

    # we want only 360 rays
    azi = azi * 360. / 2 ** azidepth
    # get azimuth of first ray
    zero_azi = int(round(azi[0]))
    # realign azimuth array to have 0deg as first ray
    azi = np.roll(azi, zero_azi, axis=0)

    # append slice elevation angle to elevs
    el = float(slices[nofslice]['posangle'])
    elevs.append(el)

    # stoprange is valid for all slices
    if nofslice == 0:
        stoprange = float(slices[nofslice]['stoprange'])
    # whereas rangestep is available in every slice
    rangestep = float(slices[nofslice]['rangestep'])

    # create range array and have it in meter
    r = np.arange(0 + rangestep/2, stoprange, rangestep) * 1000.

    # get data, make sure to only get the first 360 rays
    data_ = slices[nofslice]['slicedata']['rawdata']['data'][0:360,:]
    datadepth = float(slices[nofslice]['slicedata']['rawdata']['@depth'])
    datamin = float(slices[nofslice]['slicedata']['rawdata']['@min'])
    datamax = float(slices[nofslice]['slicedata']['rawdata']['@max'])
    data_ = datamin + data_ * (datamax - datamin) / 2 ** datadepth

    # make sure the data is aligned to zero azimuth == due north
    data_ = np.roll(data_, zero_azi, axis=0)


S = []
h = []
S2 = []
el1 = np.deg2rad(0.2)
for n in range(500):
    S_ = 0.250 + 0.500*n
    h_ = (S_**2 + ((4/3)*6371)**2 + 2*S_*(4/3)*6371*np.sin(el1))**0.5 -(4/3)*6371 + 0.053
    S2_ = (4/3)*6371*np.arcsin(S_*np.cos(el1)/((4/3)*6371+0.053))

    S.append(S_)   
    h.append(h_)
    S2.append(S2_)

S = np.array(S)
h = np.array(h)
S2 = np.array(S2)

alpha = S2/6371

lam=[]
psi=[]

for i in azi:
    lam_ = 13.072833 + (180/(np.pi))*alpha*np.cos(np.deg2rad(i))
    psi_ = 80.288333 + (180/(np.pi))*np.arcsin((np.sin(np.deg2rad(i))*np.sin(S2/6371))/(np.cos(lam_*(np.pi)/180)))
    lam.append(lam_)
    psi.append(psi_)
   
  
    lat = np.array(lam)
    lon = np.array(psi)
       

lat1 = np.reshape(lat, (1,np.product(lat.shape)))
lon1 = np.reshape(lon, (1,np.product(lon.shape)))


lat1 and lon1 are the values of coordinates for each range bins

 

[Syamsul Bahri]

Hi Krisna
what this value S_ = 0.250 + 0.500*n , 0.250 and 0.500 ?

[Kai Muehlbauer]

Hi Krishna and Syamsul,

thanks for sharing the code, Krishna. And sorry for being late, this got
somehow buried.

TL;DR: There is no need to calculate lon,lat from scratch, just use
wradlib functions. A fast basemap replacement is shown.

Since you are trying to plot on some map, the basemap package is not
supported by wradlib. We developed own functionality to make use of DEM
raster data and vector shapefiles to plot data on nice looking maps.

A recent gist for creating DEM underlay is available with this gist
notebook here:

https://gist.github.com/wradlib/8a2e10ff6c2f255256e00b8d2fb43254

How to plot data (borders, rivers, points) on the map is described in
this example section:

http://wradlib.org/wradlib-docs/latest/notebooks/visualisation/wradlib_overlay.html

To adapt/convert your radar data to longitude, latitude (e.g. wgs84) you
also can take advantage of wradlib functions. Just take your range (r)
and azimuth (az) arrays of one scan, the corresponding scan elevation
(el) and the sitecoords:

lon, lat, alt = wradlib.georef.polar2lonlatalt_n(r, az, el, sitecoords)

You'll get the lon, lat, alt for each of your scan bins. For more
parameters (earth radius, refractivity) see also:

http://wradlib.org/wradlib-docs/latest/generated/wradlib.georef.polar2lonlatalt_n.html#wradlib.georef.polar2lonlatalt_n

If you want to use the function wradlib.vis.plot_max_plan_and_vert():

http://wradlib.org/wradlib-docs/latest/generated/wradlib.vis.plot_max_plan_and_vert.html

you can have a look into this recipe:

http://wradlib.org/wradlib-docs/latest/notebooks/workflow/recipe2.html

There you would just adapt the recipe to your data and change the
projection proj.ImportFromEPSG(32632) to proj.ImportFromEPSG(4326) for
wgs84.

HTH,
Kai

Am 04.08.2016 um 07:57 schrieb Krishna chandramouli:
> Hi,
>
> I have been trying to plot the radar reflectivities read from Gematronicks
> radar .vol data.
>
> Thanks to Kai for his help succeeding it.
>
> I have referred some manual and i manually calculated the lat/lon/height
> corresponding to each reflectivity data from azimuth, elevation and range
> bins.
>
> Here is the code i used .
>
> I need to plot the data in the same way as the function
> vis.plot_max_plan_and_vert but with lat and lon ordinates
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>

> *import wradlib as wrlfrom mpl_toolkits.basemap import Basemapimport numpy
> as npimport matplotlib.pyplot as pltfrom osgeo import osrfilename =
> '/home/krishna/Downloads/IITM_DATA/rawdata/2015120123505200dBZ.vol'rbdict =
> wrl.io.read_Rainbow(filename)lon = rbdict['volume']['sensorinfo']['lon']lat

> = rbdict['volume']['sensorinfo']['lat']sitecoords = (float(lat),

> float(lon), 53.)elevs = []numele =
> int(rbdict['volume']['scan']['pargroup']['numele'])slices =

> rbdict['volume']['scan']['slice']proj =

> osr.SpatialReference()proj.ImportFromEPSG(24383)for nofslice in

> range(numele): # there are slices with 361 rays # make sure to only
> get first 360 rays azi =
> slices[nofslice]['slicedata']['rayinfo']['data'][0:360] azidepth =
> float( slices[nofslice]['slicedata']['rayinfo']['@depth']) # we
> want only 360 rays azi = azi * 360. / 2 ** azidepth # get azimuth of
> first ray zero_azi = int(round(azi[0])) # realign azimuth array to
> have 0deg as first ray azi = np.roll(azi, zero_azi, axis=0) # append
> slice elevation angle to elevs el =
> float(slices[nofslice]['posangle']) elevs.append(el) # stoprange is
> valid for all slices if nofslice == 0: stoprange =
> float(slices[nofslice]['stoprange']) # whereas rangestep is available in
> every slice rangestep = float(slices[nofslice]['rangestep']) # create
> range array and have it in meter r = np.arange(0 + rangestep/2,
> stoprange, rangestep) * 1000. # get data, make sure to only get the
> first 360 rays data_ =
> slices[nofslice]['slicedata']['rawdata']['data'][0:360,:] datadepth =
> float(slices[nofslice]['slicedata']['rawdata']['@depth']) datamin =
> float(slices[nofslice]['slicedata']['rawdata']['@min']) datamax =
> float(slices[nofslice]['slicedata']['rawdata']['@max']) data_ = datamin
> + data_ * (datamax - datamin) / 2 ** datadepth # make sure the data is
> aligned to zero azimuth == due north data_ = np.roll(data_, zero_azi,

> axis=0)S = []h = []S2 = []el1 = np.deg2rad(0.2)for n in range(500): S_ =

> 0.250 + 0.500*n h_ = (S_**2 + ((4/3)*6371)**2 +
> 2*S_*(4/3)*6371*np.sin(el1))**0.5 -(4/3)*6371 + 0.053 S2_ =
> (4/3)*6371*np.arcsin(S_*np.cos(el1)/((4/3)*6371+0.053)) S.append(S_)
> h.append(h_) S2.append(S2_)S = np.array(S) h = np.array(h)S2 =

> np.array(S2)alpha = S2/6371lam=[]psi=[]for i in azi: lam_ = 13.072833 +

> (180/(np.pi))*alpha*np.cos(np.deg2rad(i)) psi_ = 80.288333 +
> (180/(np.pi))*np.arcsin((np.sin(np.deg2rad(i))*np.sin(S2/6371))/(np.cos(lam_*(np.pi)/180)))
> lam.append(lam_) psi.append(psi_) lat = np.array(lam) lon =
> np.array(psi) lat1 = np.reshape(lat, (1,np.product(lat.shape)))lon1

> = np.reshape(lon, (1,np.product(lon.shape)))*lat1 and lon1 are the values

> of coordinates for each range bins
>
>
>

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

[Syamsul Bahri]

Tanks very much Kai,
I have try like script in below but I get error :
x = arc * np.cos(np.radians(90 - az))
ValueError: operands could not be broadcast together with shapes (480,) (327,)

script:

import wradlib as wrl

import numpy as np

import matplotlib.pyplot as pl

pl.interactive(True)

import datetime as dt

from osgeo import osr

raw = wrl.io.read_Rainbow('2016022605200600dBZ.vol')

slice=0

el=float(raw['volume']['scan']['slice'][slice]['posangle'])

sitecoords = float (raw['volume']['radarinfo']['@lon']),float (raw['volume']['radarinfo']['@lat']), float (raw['volume']['radarinfo']['@alt'])

proj = osr.SpatialReference()

proj.ImportFromEPSG(4326)#23838)

azi =raw['volume']['scan']['slice'][slice]['slicedata']['rayinfo']['data']

azidepth =float(raw['volume']['scan']['slice'][slice]['slicedata']['rayinfo']['@depth'])

azirange = float(raw['volume']['scan']['slice'][slice]['slicedata']['rayinfo']['@rays'])

az = (azi * azirange / 2**azidepth)*1.1

stoprange = float(raw['volume']['scan']['slice'][slice]['stoprange'])

rangestep = float(raw['volume']['scan']['slice'][slice]['rangestep'])*1000

r = np.arange(0,stoprange*1000,rangestep)

lon, lat, alt = wrl.georef.polar2lonlatalt_n(r, az, el, sitecoords,re=6371,ke=1.333333)

On Thursday, August 4, 2016 at 12:57:27 PM UTC+7, Krishna chandramouli wrote: