# Computing cartesian and geographical coordinates for polar data¶

[1]:

import numpy as np
import wradlib.georef as georef
import wradlib.io as io
import wradlib.util as util
import warnings

warnings.filterwarnings("ignore")

/home/runner/micromamba-root/envs/wradlib-notebooks/lib/python3.10/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm


## Read the data¶

Here, we use an OPERA hdf5 dataset.

[2]:

filename = "hdf5/20130429043000.rad.bewid.pvol.dbzh.scan1.hdf"
filename = util.get_wradlib_data_file(filename)
pvol = io.read_opera_hdf5(filename)


## Count the number of datasets¶

[3]:

ntilt = 1
for i in range(100):
try:
pvol["dataset%d/what" % ntilt]
ntilt += 1
except Exception:
ntilt -= 1
break


## Define radar location and scan geometry¶

[4]:

nrays = int(pvol["dataset1/where"]["nrays"])
nbins = int(pvol["dataset1/where"]["nbins"])
rscale = int(pvol["dataset1/where"]["rscale"])
coord = np.empty((ntilt, nrays, nbins, 3))
for t in range(ntilt):
elangle = pvol["dataset%d/where" % (t + 1)]["elangle"]
coord[t, ...] = georef.sweep_centroids(nrays, rscale, nbins, elangle)
sitecoords = (pvol["where"]["lon"], pvol["where"]["lat"], pvol["where"]["height"])
print(coord.shape)

(5, 360, 960, 3)


## Retrieve azimuthal equidistant coordinates and projection¶

[5]:

coords, proj_radar = georef.spherical_to_xyz(
coord[..., 0], coord[..., 1], coord[..., 2], sitecoords, squeeze=True
)
test = coords[0, 90, 0:960:60, 0]
print(test)

[1.24984800e+02 1.51230048e+04 3.01206529e+04 4.51178353e+04
6.01144584e+04 7.51104287e+04 9.01056525e+04 1.05100036e+05
1.20093487e+05 1.35085910e+05 1.50077213e+05 1.65067302e+05
1.80056083e+05 1.95043465e+05 2.10029352e+05 2.25013652e+05]


## Retrieve geographic coordinates (longitude and latitude)¶

### Using convenience function spherical_to_proj.¶

[6]:

lonlatalt = georef.spherical_to_proj(
coord[..., 0], coord[..., 1], coord[..., 2], sitecoords
)
test = lonlatalt[0, 90, 0:960:60, 0]
print(test)

[5.50734017 5.71615328 5.92494778 6.13371911 6.34246274 6.55117413
6.75984874 6.96848204 7.17706951 7.38560664 7.5940889  7.8025118
8.01087084 8.21916154 8.42737941 8.63552   ]


### Using reproject¶

[7]:

lonlatalt1 = georef.reproject(
coords,
projection_source=proj_radar,
projection_target=georef.get_default_projection(),
)

test = lonlatalt1[0, 90, 0:960:60, 0]
print(test)

[5.50734017 5.71615328 5.92494778 6.13371911 6.34246274 6.55117413
6.75984874 6.96848204 7.17706951 7.38560664 7.5940889  7.8025118
8.01087084 8.21916154 8.42737941 8.63552   ]