Typical workflow

[DavidQ]

Hello, I have been following the tutorial here: https://debug-docs.readthedocs.io/en/test/notebooks/basics/wradlib_workflow.html#Georeferencing-and-Projection (I know it is a test, but it seems to be what I need). But even copying the examples of the tutorial, this line [16] fails:

gridded = wrl.comp.togrid(xy, grid_xy, 128000., np.array([x.mean(), y.mean()]), data.ravel(), wrl.ipol.Nearest)

due to error in the shapes:

AssertionError: Length of value array 230400 does not correspond to number of source points 46080

What I have and have thought: a file from iris radar, SRI variable. It is a 480x480 matrix of data, supposedly azimuthal equidistant. nbins = 300 with 1 Km among them.

Is it possible to do this to paint georreferencing?: (Ideas taken from the docs, too. I create a 480x480 array of data in polar coordinates, taking into account the separation of 1000 meters):

# create range and azimuth arrays accordingly
r = np.arange(0, data.shape[1], dtype=np.float)
r += (r[1] - r[0]) / 2.
r *= 1000.
az = np.arange(0, data.shape[0], dtype=np.float)
az += (az[1] - az[0]) / 2.

polargrid = np.meshgrid(r, az)
radar_location = (344.005, 28.8744, 1517)
coords, rad = wrl.georef.spherical_to_xyz(polargrid[0], polargrid[1],
                                          elevation, radar_location)
x = coords[..., 0]
y = coords[..., 1]

projection = ccrs.AzimuthalEquidistant(central_longitude=344.8, central_latitude=28)
ax = plt.axes(projection=projection)
ax.coastlines("10m")
ax.contourf(x, y, data)
plt.show()

Would this approach be valid for my problem? Apparently it is plotting reasonable stuff, but I am not sure and have no experience on this...

[Kai Mühlbauer]

Hi David,

welcome to the wradlib community.

Answers and comments follow inline below.

Am 01.05.21 um 02:19 schrieb DavidQ:

>
> Hello, I have been following the tutorial here:
> https://debug-docs.readthedocs.io/en/test/notebooks/basics/wradlib_workflow.html#Georeferencing-and-Projection

First, please note, the wradlib documentation is available here:

https://docs.wradlib.org

The particular notebook you mention is this:

https://docs.wradlib.org/en/stable/notebooks/basics/wradlib_workflow.html#Georeferencing-and-Projection

> (I know it is a test, but it seems to be what I need). But even copying the
> examples of the tutorial, this line [16] fails:
>
> gridded = wrl.comp.togrid(xy, grid_xy, 128000., np.array([x.mean(), y.mean
> ()]), data.ravel(), wrl.ipol.Nearest)
>
> due to error in the shapes:
> AssertionError: Length of value array 230400 does not correspond to number
> of source points 46080

This might be a consequence of the outdated example. I'll need to remove
this old version from RTD. Thanks for catching this.

> *What I have and have thought:* a file from iris radar, SRI variable. It is

> a 480x480 matrix of data, supposedly azimuthal equidistant. nbins = 300
> with 1 Km among them.

I think that SRI is already a gridded product. That means the data of
shape 480x480 is in some cartesian grid in a particular projection. You
should get an immediate impression by just plotting the array with
matplotlib imshow.

> Is it possible to do this to paint georreferencing?: (Ideas taken from the
> docs, too. I create a 480x480 array of data in polar coordinates, taking
> into account the separation of 1000 meters):

[code snipped]

> Would this approach be valid for my problem? Apparently it is plotting
> reasonable stuff, but I am not sure and have no experience on this...

No, this won't work if the data is already in cartesian grid. Somewhere
in the output dictionary the x and y dimensions (probably the extent,
xmin, xmax, ymin, ymax or the like) should be available. Good chance
that more metadata is available describing a particular projection.

So it might be easier to georeference as with polar data. You could use
`wradlib.io.read_iris(filename, loaddata=False)` and show the output.
Then I could suggest a way how to plot the georeferenced data. You might
send it to me off-list, if you do not want to have the metadata in the
public.

HTH!

Cheers,
Kai

[DavidQ]

Hi and thanks, Kai.

It seems I don't have permission to write only to you to show you the metadata. There are mentions to "standard parallels 1 and 2", which makes me think in a possible? PlateCarree (regular latlon) projection. I can see ('projection_type', 0), but the name does not appear ('projection_name', ''). Also of interest are the keys ('x_location', 240000), ('y_location', 240000), in a scale ('x_scale', 100000), ('y_scale', 100000), which in my matrix of 480x480, could indicate the position of the radar.

I have also painted the data and it looks like a circle in violet with the data (precipitation, I guess) in light blue, and the values out of the circle in yellow, filling the square 480x480 matrix. I can show you the plot if you grant me permission to write specifically to you.

Best, David Q.

[DavidQ]

This is the solution (below) to the problem, thanks to the immense help and support of Kai. First thing to notice is that SRI was in Azimuthal Equidistant, as it was found examining the metadata (wrl.io.read_iris(filename, loaddata=False) and using the IRIS information ( ftp://ftp.sigmet.com/outgoing/manuals/IRIS_Programmers_Manual.pdf)

import cartopy.crs as ccrs
import matplotlib.pyplot as plt

import wradlib as wrl

import os

os.environ["WRADLIB__DATA"] = "."

iris0 = wrl.io.read_iris(filename)

x_size = iris0['product_hdr']["product_configuration"]["x_size"]
y_size = iris0['product_hdr']["product_configuration"]["y_size"]

x_scale = iris0['product_hdr']["product_configuration"]["x_scale"] /100. #note the units in metadata
y_scale = iris0['product_hdr']["product_configuration"]["y_scale"] /100.
x_location = iris0['product_hdr']["product_configuration"]["x_location"]// 1000 #again, for scale reasons
y_location = iris0['product_hdr']["product_configuration"]["y_location"]// 1000

lat0 = iris0['product_hdr']["product_end"]["latitude_projection"]
lon0 = iris0['product_hdr']["product_end"]["longitude_projection"]
globe = ccrs.Globe(semimajor_axis = 6378137, inverse_flattening = 298.257224) #I needed a globe since I wasn't exactly in WGS-84. semimajor_axis has been #taken from radius of the Earth and convert it to meters
map_proj = ccrs.AzimuthalEquidistant(central_latitude=lat0, central_longitude=lon0, globe=globe)

x_dim = (np.arange(x_size) - x_location) * x_scale
y_dim = (np.arange(y_size) - y_location) * y_scale

fig = plt.figure(figsize=(20,16))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree(globe=globe))
ax.set_aspect('equal')

X, Y = np.meshgrid(x_dim, y_dim)

pm = ax.pcolormesh(X, Y, iris0['data'][0][0], transform=map_proj)
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.colorbar(pm)

plt.show()

[Kai Muehlbauer]

Hi David,

thanks for the example and the explanations.

I'll only add that users might need to use a ftp client with user
"anonymous" (and maybe also password "anonymous") to gain access to the
FTP. Below is the name of the programmers guide as of now.

ftp://ftp.sigmet.com/outgoing/manuals/IRIS-Programming-Guide-M211318EN.pdf

Cheers,
Kai

Am 05.05.21 um 10:51 schrieb DavidQ:

> This is the solution (below) to the problem, thanks to the immense help
> and support of Kai. First thing to notice is that SRI was in Azimuthal
> Equidistant, as it was found examining the metadata
> (wrl.io.read_iris(filename, loaddata=False) and using the IRIS
> information (
> ftp://ftp.sigmet.com/outgoing/manuals/IRIS_Programmers_Manual.pdf

> <ftp://ftp.sigmet.com/outgoing/manuals/IRIS_Programmers_Manual.pdf>)

> https://docs.wradlib.org <https://docs.wradlib.org>

> --
> 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>.
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/wradlib-users/99bfd1da-a60e-4d54-9fcc-805ee9da72afn%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/99bfd1da-a60e-4d54-9fcc-805ee9da72afn%40googlegroups.com?utm_medium=email&utm_source=footer>.

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