ERRORS ON GRIDDING DATA - UNEXPECTED VALUES BEYOND THE RADAR RANGE
46 views
Skip to first unread message
Andrea Antonini
Jun 7, 2021, 10:07:24 AM6/7/21
to wradlib-users
Hello everybody,
I'm new in this group. I'm trying to implement some parts of the dataflow described in:
to convert some radar polar volumes in a geotiff image (lat lon points).
All the steps seem to work correctly, but there are some problems when the data are gridded in the final lat lon grid using the functions:
wrl.comp.togrid
np.ma.masked_invalid
Some problems occurs when at the very end of the beam, the last cell has a reflectivity value different from the noise level. The interpolation process of the wrl.comp.togrid function can cause an error introducing some reflectivity values beyond the radar range, different from the nodata value.
I put in attachment two images: the first (Figure_A) is a PPI image (without gridding processing) and the second (Figure_B) after the gridding processing. You can see in the bottom left part of the image an artificial conical echo reflectivity.
Do you have any suggestion about how to correct this?
Thank you very much
Best regards
Andrea Antonini
Kai Muehlbauer
Jun 7, 2021, 10:21:12 AM6/7/21
to wradli...@googlegroups.com
Hi Andrea,
welcome to wradlib.
Can you just show the complete codelines? It looks like the data of the
last range grid is extrapolated (nearest) outwards. I suspect you create
a cartesian grid which is slightly larger than the ppi max range? Could
you lower the radius value (in togrid) accordingly?
HTH,
Kai
Am 07.06.21 um 16:07 schrieb Andrea Antonini:
> --
> 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/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com?utm_medium=email&utm_source=footer>.
welcome to wradlib.
Can you just show the complete codelines? It looks like the data of the
last range grid is extrapolated (nearest) outwards. I suspect you create
a cartesian grid which is slightly larger than the ppi max range? Could
you lower the radius value (in togrid) accordingly?
HTH,
Kai
Am 07.06.21 um 16:07 schrieb Andrea Antonini:
> 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/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com?utm_medium=email&utm_source=footer>.
Andrea Antonini
Jun 7, 2021, 10:28:31 AM6/7/21
to wradli...@googlegroups.com
Dear Kai,
thenk you for your e-mail,
following in attachment the python code that take data and plot the two images.
# Thake offset and gain info from file
gain=fcontent['dataset%d/data1/what' %i]['gain']
offset=fcontent['dataset%d/data1/what' %i]['offset']
#Take data from file
data = fcontent['dataset%d/data1/data' %i]
data = gain*np.array(data)+offset
fig = pl.figure(figsize=(24,24), dpi=100)
ax, pm = wrl.vis.plot_ppi(data, fig=fig)
pl.title('Reflectivity Map')
pl.show()
'''
#Clutter statistics
clutter = wrl.clutter.filter_gabella(data, tr1=12, n_p=6, tr2=1.1)
pl.figure(figsize=(10,10))
ax, pm = wrl.vis.plot_ppi(clutter, cmap=pl.cm.gray)
pl.title('Clutter Map')
#Decluttering
data_no_clutter = wrl.ipol.interpolate_polar(data, clutter)
pl.figure(figsize=(10,10))
ax, pm = wrl.vis.plot_ppi(data_no_clutter) # simple diagnostic plot
cbar = pl.colorbar(pm, shrink=0.75)
'''
elev = "elev%d" % i
dataset = "dataset%d" % i
shapefile_name = "%s_centroidi%d" % (hdf5_file_base , i)
#print(elev)
# Single elevation geotiff file name
geotiff_elev_name = "%s_%s.tif" % (hdf5_file_base , elev)
#print(geotiff_elev_name)
geotiff_list.append("%s%s" % (wdir,geotiff_elev_name))
#Get the elevation angle
elangle = fcontent['dataset%d/where' % i]['elangle']
#print(elangle)
'''
coords, proj_radar = georef.spherical_to_xyz(
coord[..., 0],
np.degrees(coord[..., 1]),
coord[..., 2],
sitecoords
)
'''
#Manipulating coordinates
# Get the x,y,z coordinatres of each point of the scan. coords dimension is numelev x numazimuth x numnbeams x 3
coords, proj_radar = wrl.georef.spherical_to_xyz(polargrid[0], polargrid[1],
elangle, sitecoords)
#Get the final grid coordinates
#lonlatalt coordinates are in lat lon height lonlatalt dimension is numazimuth x numnbeams x 3
lonlatalt = georef.spherical_to_proj(polargrid[0], polargrid[1],
elangle, sitecoords)
#lonlatalt1 dimension is is numelev x numazimuth x numnbeams x 3
lonlatalt1 = georef.reproject(
coords,
projection_source=proj_radar,
projection_target=georef.get_default_projection()
)
#Coordinates management
x = lonlatalt1[..., 0]
y = lonlatalt1[..., 1]
#cent_lon will be used in exporting the geotiff
cent_lon = lonlatalt[..., 0:2]
#Setup the grids
xgrid = np.linspace(x.min(), x.max(), 2400)
ygrid = np.linspace(y.min(), y.max(), 2400)
grid_xy = np.meshgrid(xgrid, ygrid)
grid_xy = np.vstack((grid_xy[0].ravel(), grid_xy[1].ravel())).transpose()
x_min = x.min()
y_max = y.max()
x_max = x.max()
y_min = y.min()
x_ravel = x.ravel()
x_ravel_none = x.ravel()[:,None]
y_ravel = y.ravel()
y_ravel_none = y.ravel()[:,None]
xy=np.concatenate([x.ravel()[:,None],y.ravel()[:,None]], axis=1)
print(xy.shape)
print(x_ravel_none.shape)
print(y_ravel_none.shape)
print(wrl.__version__)
# Now we transfer the polar data to the grid
gridded = wrl.comp.togrid(xy, grid_xy, 108000., np.array([x.mean(), y.mean()]), data.ravel(), wrl.ipol.Nearest)
#gridded = np.ma.masked_where(gridded == 0, gridded).reshape((len(xgrid), len(ygrid)))
gridded = np.ma.masked_invalid(gridded).reshape((len(xgrid), len(ygrid)))
#Printing the geolocated data
fig = pl.figure(figsize=(10,8))
ax = pl.subplot(111, aspect="equal")
pm = pl.pcolormesh(xgrid, ygrid, gridded)
pl.colorbar(pm, shrink=0.75)
pl.xlabel("Easting (m)")
pl.ylabel("Northing (m)")
pl.xlim(min(xgrid), max(xgrid))
pl.ylim(min(ygrid), max(ygrid))
pl.show()
gain=fcontent['dataset%d/data1/what' %i]['gain']
offset=fcontent['dataset%d/data1/what' %i]['offset']
#Take data from file
data = fcontent['dataset%d/data1/data' %i]
data = gain*np.array(data)+offset
fig = pl.figure(figsize=(24,24), dpi=100)
ax, pm = wrl.vis.plot_ppi(data, fig=fig)
pl.title('Reflectivity Map')
pl.show()
'''
#Clutter statistics
clutter = wrl.clutter.filter_gabella(data, tr1=12, n_p=6, tr2=1.1)
pl.figure(figsize=(10,10))
ax, pm = wrl.vis.plot_ppi(clutter, cmap=pl.cm.gray)
pl.title('Clutter Map')
#Decluttering
data_no_clutter = wrl.ipol.interpolate_polar(data, clutter)
pl.figure(figsize=(10,10))
ax, pm = wrl.vis.plot_ppi(data_no_clutter) # simple diagnostic plot
cbar = pl.colorbar(pm, shrink=0.75)
'''
elev = "elev%d" % i
dataset = "dataset%d" % i
shapefile_name = "%s_centroidi%d" % (hdf5_file_base , i)
#print(elev)
# Single elevation geotiff file name
geotiff_elev_name = "%s_%s.tif" % (hdf5_file_base , elev)
#print(geotiff_elev_name)
geotiff_list.append("%s%s" % (wdir,geotiff_elev_name))
#Get the elevation angle
elangle = fcontent['dataset%d/where' % i]['elangle']
#print(elangle)
'''
coords, proj_radar = georef.spherical_to_xyz(
coord[..., 0],
np.degrees(coord[..., 1]),
coord[..., 2],
sitecoords
)
'''
#Manipulating coordinates
# Get the x,y,z coordinatres of each point of the scan. coords dimension is numelev x numazimuth x numnbeams x 3
coords, proj_radar = wrl.georef.spherical_to_xyz(polargrid[0], polargrid[1],
elangle, sitecoords)
#Get the final grid coordinates
#lonlatalt coordinates are in lat lon height lonlatalt dimension is numazimuth x numnbeams x 3
lonlatalt = georef.spherical_to_proj(polargrid[0], polargrid[1],
elangle, sitecoords)
#lonlatalt1 dimension is is numelev x numazimuth x numnbeams x 3
lonlatalt1 = georef.reproject(
coords,
projection_source=proj_radar,
projection_target=georef.get_default_projection()
)
#Coordinates management
x = lonlatalt1[..., 0]
y = lonlatalt1[..., 1]
#cent_lon will be used in exporting the geotiff
cent_lon = lonlatalt[..., 0:2]
#Setup the grids
xgrid = np.linspace(x.min(), x.max(), 2400)
ygrid = np.linspace(y.min(), y.max(), 2400)
grid_xy = np.meshgrid(xgrid, ygrid)
grid_xy = np.vstack((grid_xy[0].ravel(), grid_xy[1].ravel())).transpose()
x_min = x.min()
y_max = y.max()
x_max = x.max()
y_min = y.min()
x_ravel = x.ravel()
x_ravel_none = x.ravel()[:,None]
y_ravel = y.ravel()
y_ravel_none = y.ravel()[:,None]
xy=np.concatenate([x.ravel()[:,None],y.ravel()[:,None]], axis=1)
print(xy.shape)
print(x_ravel_none.shape)
print(y_ravel_none.shape)
print(wrl.__version__)
# Now we transfer the polar data to the grid
gridded = wrl.comp.togrid(xy, grid_xy, 108000., np.array([x.mean(), y.mean()]), data.ravel(), wrl.ipol.Nearest)
#gridded = np.ma.masked_where(gridded == 0, gridded).reshape((len(xgrid), len(ygrid)))
gridded = np.ma.masked_invalid(gridded).reshape((len(xgrid), len(ygrid)))
#Printing the geolocated data
fig = pl.figure(figsize=(10,8))
ax = pl.subplot(111, aspect="equal")
pm = pl.pcolormesh(xgrid, ygrid, gridded)
pl.colorbar(pm, shrink=0.75)
pl.xlabel("Easting (m)")
pl.ylabel("Northing (m)")
pl.xlim(min(xgrid), max(xgrid))
pl.ylim(min(ygrid), max(ygrid))
pl.show()
You received this message because you are subscribed to a topic in the Google Groups "wradlib-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/wradlib-users/4LXXzFAUPWs/unsubscribe.
To unsubscribe from this group and all its topics, send an email to wradlib-user...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/wradlib-users/727aa95c-9234-a023-cab8-59713a247b6a%40uni-bonn.de.
Kai Muehlbauer
Jun 7, 2021, 10:59:32 AM6/7/21
to wradli...@googlegroups.com
Dear Andrea,
my first guess would be to lower the "108000." (in the call to togrid)
gradually until you see the erroneous extrapolation reduced as wanted.
But this might have side effects.
Another possibility is to not use latlon but some other projection (eg.
azimuthal equidistant, UTM) to create the grid and reproject afterwards
to latlon.
HTH,
Kai
Am 07.06.21 um 16:28 schrieb Andrea Antonini:
> <kai.mue...@uni-bonn.de <mailto:kai.mue...@uni-bonn.de>> ha scritto:
> <mailto:wradlib-users%2Bunsu...@googlegroups.com>
> > <mailto:wradlib-user...@googlegroups.com
> <mailto:wradlib-users%2Bunsu...@googlegroups.com>>.
> <https://groups.google.com/d/msgid/wradlib-users/727aa95c-9234-a023-cab8-59713a247b6a%40uni-bonn.de>.
>
> --
> <https://groups.google.com/d/msgid/wradlib-users/CABHZe%2B5-VJvxbvnq7bFxwvgB54usWMTRviFHGL_vFej8v3WcFA%40mail.gmail.com?utm_medium=email&utm_source=footer>.
my first guess would be to lower the "108000." (in the call to togrid)
gradually until you see the erroneous extrapolation reduced as wanted.
But this might have side effects.
Another possibility is to not use latlon but some other projection (eg.
azimuthal equidistant, UTM) to create the grid and reproject afterwards
to latlon.
HTH,
Kai
Am 07.06.21 um 16:28 schrieb Andrea Antonini:
> > <mailto:wradlib-user...@googlegroups.com
> <mailto:wradlib-users%2Bunsu...@googlegroups.com>>.
> > To view this discussion on the web, visit
> >
> https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com>
>
> >
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com?utm_medium=email&utm_source=footer
> >
> https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com>
>
> >
> <https://groups.google.com/d/msgid/wradlib-users/5b462c5c-4674-4fa1-bbc0-e6f9b3a3d815n%40googlegroups.com?utm_medium=email&utm_source=footer>>.
>
> --
> You received this message because you are subscribed to a topic in
> the Google Groups "wradlib-users" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/wradlib-users/4LXXzFAUPWs/unsubscribe
> <https://groups.google.com/d/topic/wradlib-users/4LXXzFAUPWs/unsubscribe>.
>
> --
> You received this message because you are subscribed to a topic in
> the Google Groups "wradlib-users" group.
> To unsubscribe from this topic, visit
> https://groups.google.com/d/topic/wradlib-users/4LXXzFAUPWs/unsubscribe
> To unsubscribe from this group and all its topics, send an email to
> wradlib-user...@googlegroups.com
> <mailto:wradlib-users%2Bunsu...@googlegroups.com>.
> wradlib-user...@googlegroups.com
> To view this discussion on the web, visit
> https://groups.google.com/d/msgid/wradlib-users/727aa95c-9234-a023-cab8-59713a247b6a%40uni-bonn.de
> <https://groups.google.com/d/msgid/wradlib-users/727aa95c-9234-a023-cab8-59713a247b6a%40uni-bonn.de>.
>
> --
> 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/CABHZe%2B5-VJvxbvnq7bFxwvgB54usWMTRviFHGL_vFej8v3WcFA%40mail.gmail.com
> 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/CABHZe%2B5-VJvxbvnq7bFxwvgB54usWMTRviFHGL_vFej8v3WcFA%40mail.gmail.com?utm_medium=email&utm_source=footer>.
Reply all
Reply to author
Forward
0 new messages