DX data to lat/lon regular grid

[Franck Chopin]

Good morning

 

I am currently working (and this is completely new to me) with DX

radar data from DWD.

 

In order to read this data I have installed the wradlib python library

into a Ubuntu workstation.

 

The good thing with it is that I am now learning a programing language

which is completely new to me.

 

The bad thing with it is that I am struggling a lot !

 

My intent is to read the DX data, to convert it into a regular  grid

where each grid cell is 0.005 degrees x 0.005 degrees and to export

this data into a netcdf file which is easier to manipulate for me.

 

From various websites I looked at, it seems that doing this data

manipulation is feasible via the wradlib library.

 

The code I am currently writing is the following coming mainly

from the http://wradlib.bitbucket.org/ website.

 

*******

import wradlib

import pylab as pl
data, metadata = wradlib.io.readDX("10908DX130728_1635.ARC")
import numpy as np
radar_location = (8.005, 47.8744, 1517) # (lon, lat, alt) in decimal degree and meters # FBG Radar
elevation = 0.5 # in degree
azimuths = np.arange(0,360) # in degrees
ranges = np.arange(0, 128000., 1000.) # in meters
polargrid = np.meshgrid(ranges, azimuths)
lon, lat, alt = wradlib.georef.polar2lonlatalt_n(polargrid[0], polargrid[1], elevation, radar_location)
# here lon, lat and alt are 3 arrays of 360 by 128

# At this stage we have for each bin of the radar its longitude, its latitude and its altitude
#The objective is to regrid it to a regular grid as needed: a regular grid (latitude/longitude) where each grid cell is 0.005 degrees x 0.005 degrees.

longrid = np.linspace(lon.min(), lon.max(), (lon.max()-lon.min())/0.005)
latgrid = np.linspace(lat.min(), lat.max(), (lat.max()-lat.min())/0.005)
# creation of 2 nd arrays of dimension 1 of lineraly spaced lon values et lineraly spaced lat values with an offset of 0.005 degrees

grid_lonlat = np.meshgrid(longrid, latgrid)

grid_lonlat = np.vstack((grid_lonlat[0].ravel(), grid_lonlat[1].ravel())).transpose()

xy=np.concatenate([lon.ravel()[:,None],lat.ravel()[:,None]], axis=1)

#gridded = wradlib.comp.togrid(xy, grid_lonlat, 128000., np.array([lon.mean(), lat.mean()]), data.ravel(), wradlib.ipol.Nearest)

gridded = wradlib.comp.togrid(xy, grid_lonlat, 1, np.array([lon.mean(), lat.mean()]), data.ravel(), wradlib.ipol.Nearest)

gridded = np.ma.masked_invalid(gridded)
gridded=gridded.reshape((len(latgrid), len(longrid)))

fig = pl.figure(figsize=(10,8))
ax = pl.subplot(111, aspect="equal")
pm = pl.pcolormesh(longrid, latgrid, gridded)
pl.colorbar(pm, shrink=0.75)
pl.xlabel("Longitude")
pl.ylabel("Latitude")
pl.show()

*********

 

 

At this stage, I obtain the regular grid I want BUT I have strong suspicions on the projection. I am not sure I am using the wradlib.comp.togrid function the way I should (in red above). The third parameter is supposed to be the radius of the radar which is tricky. My grid

being in Lat/Lon the radar image is not a circle. I do not know if the way I am doing the projection to a regular grid is correct...

Could you please help me writing a correct code ?

 

I will next have to export the regular grid into a netcdf .grd file but before thinking about this next step I would like to be sure what I am doing now is correct.

 

Thank you for your help and best regards,

 

Franck

[Maik Heistermann]

Dear Franck,

thanks for your request. On first look, your application is appears to point in the right direction. The issue, I guess, is with your target grid definition. A lat/lon grid cannot be interpreted as a cartesian reference system, meaning that the grid spacing (in terms of real distances) is not homogeneous over your spatial domain. The function wradlib.comp.togrid, however, makes sense only in case the target coordinates represent a cartesian reference system. Otherwise, the "radius" parameter represents different distances in north and south directions, although it is intended to represent the radar range which is, of course, the same in any direction around the radar.

So, my question would be: Do you actually need a regular lat/lon grid? Or would it also help you to first grid your data on a cartesian grid, and then reproject the grid points to lat/lon? Otherwise, if you need a regular lat/long grid, I would suggest you do not use the wradlib.comp.togrid function, but to interpolate "directly" (using e.g. wradlib.ipol.interpolate) to your target lat/longrid. In that case, however, you will find interpolated values on any grid point, including those beyond the radar range. Then you would need to take case of masking these values by youself.

Does that make sense? I will try to look into the code in more detail - maybe I can come up with a tangible solution.

Bests regards,
Maik

[Maik Heistermann]

Dear Franck,

just as an example: This would be the code in case you would use a target grid in a cartesian reference system (here: Gauss Krueger Zone 3). Here I used DX data from radar Feldberg in SW-Germany (also contained in the wradlib "examples/data" directory.

import wradlib
import numpy as np
import pylab as pl

data, meta = wradlib.io.readDX("raa00-dx_10908-0806021655-fbg---bin")

# set scan geometry and radar coordinates
r  = np.arange(500.,128500.,1000.)
az  = np.arange(0,360)
sitecoords  = (8.005, 47.8744)

# create osr projection using epsg number for GK Zone 3
proj_gk3 = wradlib.georef.epsg_to_osr(31467)
# lat/lon coordiantes of bin centroids
cent_lon, cent_lat = wradlib.georef.polar2centroids(r, az, sitecoords)
x, y = wradlib.georef.reproject(cent_lon, cent_lat, projection_target=proj_gk3)
coords = np.array([x.ravel(),y.ravel()]).transpose()

# define target grid for composition
xmin, xmax, ymin, ymax = x.min(), x.max(), y.min(), y.max()
trgx = np.linspace(x.min(),x.max()+1000.,1000.)
trgy = np.linspace(y.min(),y.max()+1000.,1000.)
trgcoords = wradlib.util.gridaspoints(trgy,trgx)

# interpolate polar radar-data and quality data to the grid
gridded = wradlib.comp.togrid(coords, trgcoords, r.max()+500., coords.mean(axis=0), data.ravel(), wradlib.ipol.Nearest)

# Example map

fig = pl.figure(figsize=(10,8))
ax = pl.subplot(111, aspect="equal")

pm = pl.pcolormesh(trgx, trgy, np.ma.masked_invalid( gridded.reshape((len(trgx),len(trgy))) ), cmap="spectral")
grd = pl.grid()
cb = pl.colorbar(pm, shrink=0.5)
cb.set_label("dBZ")

Bests,

Maik

Am Montag, 28. September 2015 10:29:16 UTC+2 schrieb Franck Chopin:

[Maik Heistermann]

And this would be the same example, showing both the grid in Gaus Krueger Zone 3, and the same grid projected to lat/lon. However, you need to be aware that in lat/lon, this is not a regular grid anymore, but rather a point cloud (which only matplotlibs pcolormesh function helps us to represent as a grid).

import wradlib

import numpy as np

import pylab as pl

data, meta = wradlib.io.readDX("raa00-dx_10908-0806021655-fbg---bin")

# set scan geometry and radar coordinates

r  = np.arange(500.,128500.,1000.)

az  = np.arange(0,360)

sitecoords  = (8.005, 47.8744)

# create osr projection using epsg number for GK Zone 3

proj_gk3 = wradlib.georef.epsg_to_osr(31467)

# lat/lon coordiantes of bin centroids

cent_lon, cent_lat = wradlib.georef.polar2centroids(r, az, sitecoords)

x, y = wradlib.georef.reproject(cent_lon, cent_lat, projection_target=proj_gk3)

coords = np.array([x.ravel(),y.ravel()]).transpose()

# define target grid for composition

xmin, xmax, ymin, ymax = x.min(), x.max(), y.min(), y.max()

trgx = np.linspace(x.min(),x.max()+1000.,1000.)

trgy = np.linspace(y.min(),y.max()+1000.,1000.)

trgcoords = wradlib.util.gridaspoints(trgy,trgx)

# interpolate polar radar-data and quality data to the grid

gridded = wradlib.comp.togrid(coords, trgcoords, r.max()+500., coords.mean(axis=0), data.ravel(), wradlib.ipol.Nearest)

# Example map

fig = pl.figure(figsize=(10,8))

ax = pl.subplot(111, aspect="equal")

pm = pl.pcolormesh(trgx, trgy, np.ma.masked_invalid( gridded.reshape((len(trgx),len(trgy))) ), cmap="spectral")

grd = pl.grid()

cb = pl.colorbar(pm, shrink=0.5)

cb.set_label("dBZ")

trglon, trglat = wradlib.georef.reproject(trgx, trgy,

                                            projection_source=proj_gk3,

                                            projection_target=wradlib.georef.get_default_projection() )

# Example map

fig = pl.figure(figsize=(12,8))

ax = pl.subplot(121, aspect="equal")

pm = pl.pcolormesh(trgx, trgy, np.ma.masked_invalid( gridded.reshape((len(trgx),len(trgy))) ), cmap="spectral")

grd = pl.grid()

cb = pl.colorbar(pm, shrink=0.5)

cb.set_label("dBZ")

pl.title("Gauss Krueger Zone 3")

ax = pl.subplot(122, aspect="equal")

pm = pl.pcolormesh(trglon, trglat, np.ma.masked_invalid( gridded.reshape((len(trglon),len(trglat))) ), cmap="spectral")

grd = pl.grid()

cb = pl.colorbar(pm, shrink=0.5)

cb.set_label("dBZ")

pl.title("Reprojected to lat/lon")

Bests,
Maik

Am Montag, 28. September 2015 10:29:16 UTC+2 schrieb Franck Chopin:

[Franck Chopin]

Thank you very much Maik for this quick and very detailled answer. I will have a look at it to digest it.

Yes I was aware that in the lat/Lon regular grid the distances were not kept. I just hope that there were a way to set up the parameters so that it was working for lat/lons regular grids.

Thank you very much for your help ! I really apprecaite it!

Best regards,

Franck

[Franck Chopin]

Dear Maik,

Unfortunately I do need to have my final grid in a lat/lon regular grid of 0.005 degrees. This is something I cannot escape.

Maybe I can use your second code to obtain a non regular lat lon grid and then reinterpolate it. This way I might not need to use a mask.

Thank you,

Franck

Le lundi 28 septembre 2015 10:29:16 UTC+2, Franck Chopin a écrit :

[Franck Chopin]

Maik,

I have a strange answer when using your code about one attribute:

>>> proj_gk3 = wradlib.georef.epsg_to_osr(31467)

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

AttributeError: 'module' object has no attribute 'epsg_to_osr'

When looking at 

>>> dir(wradlib.georef)

this is true that the epsg_to_osr does not appear in the list:

['__builtins__', '__doc__', '__file__', '__name__', '__package__', '__pol2lonlat', '_check_polar_coords', '_doctest_', '_get_azimuth_resolution', '_get_range_resolution', '_is_sorted', '_latscale', '_lonscale', 'aeq2hor', 'apichange_kwarg', 'arc_distance_n', 'beam_height_n', 'centroid2polyvert', 'create_projstr', 'exit', 'gdal', 'get_default_projection', 'get_earth_radius', 'get_radolan_coords', 'get_radolan_grid', 'hor2aeq', 'np', 'osr', 'pixel_coordinates', 'pixel_to_map', 'pixel_to_map3d', 'polar2centroids', 'polar2lonlat', 'polar2lonlatalt', 'polar2lonlatalt_n', 'polar2polyvert', 'proj4_to_osr', 'projected_bincoords_from_radarspecs', 'read_gdal_coordinates', 'read_gdal_projection', 'read_gdal_values', 'reproject', 'sweep_centroids', 'warnings']

Now when checking the wradlib I have :

~$ pip list | grep wradlib

wradlib (0.5.1)

I do have the 0.5.1 version which I think is the latest...

If you have any suggestion, I take it !

Thank you...

Franck

Le lundi 28 septembre 2015 10:29:16 UTC+2, Franck Chopin a écrit :

[Maik Heistermann]

Dear Franck,

I see. Still, I would suggest to interpolate directly from polar (bin) coordinates to your lat/lon grid, and mask values outside the radar circle afterwards. This is computationally more efficient (otherwise, you'd have to to an interpolation twice). How about this code:

import wradlib

import numpy as np

import pylab as pl

data, meta = wradlib.io.readDX("raa00-dx_10908-0806021655-fbg---bin")

# set scan geometry and radar coordinates

r  = np.arange(500.,128500.,1000.)

az  = np.arange(0,360)

sitecoords  = (8.005, 47.8744)

# lat/lon coordiantes of bin centroids

lon, lat = wradlib.georef.polar2centroids(r, az, sitecoords)

coords = np.array([lon.ravel(),lat.ravel()]).transpose()

# define target grid (fine tune to your needs)

trglon = np.arange(lon.min(),lon.max(),0.005)

trglat = np.arange(lat.min(),lat.max(),0.005)

trgcoords = wradlib.util.gridaspoints(trglat,trglon)

# Gridding (even beyond the maximum radar range)

gridded = wradlib.ipol.interpolate(coords, trgcoords, data.ravel(), wradlib.ipol.Nearest)

# Create mask in order to "remove" values outside the radar circle

# We do this in a Cartesian reference system, e.g. GK Zone 3

proj_gk3 = wradlib.georef.epsg_to_osr(31467)

sitecoords_repr = wradlib.georef.reproject(sitecoords, projection_target=proj_gk3)

trgcoords_repr = wradlib.georef.reproject(trgcoords, projection_target=proj_gk3)

mask = np.where( ((trgcoords_repr-sitecoords_repr)**2).sum(axis=-1) > r.max()**2 )[0]

# Mask out values outside the circle

gridded[mask] = np.nan

# Example map

fig = pl.figure(figsize=(10,8))

ax = pl.subplot(111, aspect="equal")

pm = pl.pcolormesh(trglon, trglat,

                   np.ma.masked_invalid( gridded.reshape((len(trglat),len(trglon))) ),

                   cmap="spectral")

grd = pl.grid()

cb = pl.colorbar(pm, shrink=0.5)

cb.set_label("dBZ")

Bests,
Maik

[Maik Heistermann]

Yes, it might be true that this function is available only in the bleeding edge version of wradlib. I'd suggest you install that version from our onlie repository:

Just download "http://bitbucket.org/wradlib/wradlib/get/default.zip", extract, and run "python setup.py install".

Bests,
Maik

[Franck Chopin]

Thank you very much Maik!

It seems to work. I will now work on the next phase of the work. Will come back to you if I struggle but I will try to do it myself first. If not I will not learn...

Your help is much appreciated !

Best regards,

Franck