Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
lon0 = 10. # central meridian of projection
lat0 = 60. # standard parallel of projection
sinlat0 = np.sin(np.radians(lat0))
fac = (6370.040 ** 2.) * ((1. + sinlat0) ** 2.)
lon = np.degrees(np.arctan((-x / y))) + lon0
lat = np.degrees(np.arcsin((fac - (x ** 2. + y ** 2.)) /
(fac + (x ** 2. + y ** 2.))))
radolan_grid = np.dstack((lon, lat))
else:
# create radolan projection osr object
proj_stereo = projection.create_osr("dwd-radolan")
# create wgs84 projection osr object
proj_wgs = projection.get_default_projection()
radolan_grid = projection.reproject(radolan_grid,
projection_source=proj_stereo,
projection_target=proj_wgs)
return radolan_grid
# as described in the format description
phi_0 = np.radians(60)
phi_m = np.radians(lat)
lam_0 = 10
lam_m = lon
lam = np.radians(lam_m - lam_0)
er = 6370.040
m_phi = (1 + np.sin(phi_0)) / (1 + np.sin(phi_m))
x = er * m_phi * np.cos(phi_m) * np.sin(lam)
y = - er * m_phi * np.cos(phi_m) * np.cos(lam)
else:
# create radolan projection osr object
proj_stereo = projection.create_osr("dwd-radolan")
# create wgs84 projection osr object
proj_wgs = projection.get_default_projection()
x, y = projection.reproject(lon, lat, projection_source=proj_wgs,
projection_target=proj_stereo)
return x, y
9.0000, 48.0000, 0.0000
9.0000, 48.0000, 0.0000
9.0000, 48.9981, 725.7160
10.4872, 47.9904, 725.7160
9.0000, 47.0017, 725.7160
7.5131, 47.9904, 1694.2234
Here, the coordinates of the east and west directions won't come to lie on
the latitude of the site because the beam doesn't travel along the latitude
circle but along a great circle.
See :ref:`/notebooks/basics/wradlib_workflow.ipynb#\
Georeferencing-and-Projection`.
"""
if proj is None:
proj = projection.get_default_projection()
xyz, rad = spherical_to_xyz(r, phi, theta, sitecoords, re=re, ke=ke,
squeeze=True)
# reproject aeqd to destination projection
coords = projection.reproject(xyz, projection_source=rad,
projection_target=proj)
return coords