Conversion from Lat/Lon to Coordinates

Python
1

Hi,

I’m trying to determine the location of radar bins relative to a point, but I’m stuck on step [5] of this tutorial:
tutorial

sitecoords = (9.7839, 48.5861)

Those values I assume to be decimal latitude and longitude, is it right?

Then they become

#Coordinates of the rain gages in Gauss-Krueger 3 coordinates
x, y = np.array([3557880, 3557890]), np.array([5383379, 5383375])

I have the Lat, Lon values of my example point:

[-9.6616, -35.6972]

How do I get to those new coordinate values?

Thanks in advance
Dan

1 Like
2

@d-siden Two issues here:

  • wradlib uses decimal (lon, lat) (OAMS_TRADITIONAL_GIS_ORDER) not (lat,lon), you would need to specify your coordinates with (lon, lat)
  • the coordinates of the raingauges are given arbitrary and they fall into the radar domain of the specified radar

It looks like this example needs some refresh, too.

If you want to work along that example with your own radar location, you would also want to use another projection. Gauss-Krueger 3 is a projection for the region of mid-Europe (Germany).

Please ask back, if need be.

3

I’ve been having some troubles finding the nearest bins to a place in town. I think I’ve made some progress. If I don’t manage to finish this today I’ll consider bothering you again. :melting_face:

4

Now I’m curious if I got this right, because in the above mentioned tutorial, I suppose this is the amount of bins in each azimuth.

r = np.arange(1, 129)

But in step 2 of this tutorial it is different.

r = np.arange(1, 129) * 1000

My data has 1000 bins in each azimuth. I’ll try plotting with shapefiles on the background to see if I can find the perfect spot for me.

1 Like
5

My assumption that the 1000 value above is a metric range resolution seems to be correct.
In my case, my range step (range resolution) is 250 meters. Alongside with the number of bins (1000), the r argument for PolarNeighbours() function is

r=numpy.arange(1,nbins+1)*range_resolution

This resulted in a correct alignment with a shapefile in the background.

6

@d-siden Great to see you have found your way through those examples. And also good to see how the examples could be improved. If you are up to do some revision/updating I’m more than happy to help.

1 Like
7

If you want to update the tutorial, here is how it is working:

Extract bin values from a polar radar data set at rain gage locations

Suppose you already have the following radar data:

nbins: number of bins in each azimuthal angle;
range_resolution (a.k.a.: range step): distance in radius between bins (in meters);
radar_Lon, radar_Lat: radar site longitude and latitude;
azi: azimuth angles in a scan (degrees);

The azimuth angles are a numpy.ndarray object cointaining the angles in the order the file provides them; like the way it results from this tutorial.

The projection used in this case is WGS84.
This projection allows to input the point location straight as decimal degrees. Although it might not be the more appropriate.

projec = wradlib.georef.epsg_to_osr(4326)

If needed, point location can be inserted at a different projection through this process:

  1. Use the following pyproj library function to create the desired “transformer” (in this case a PROJ string is used to create WGS84, check this doc, there are many ways to do it):
import pyproj
from pyproj import Proj
transf = Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
  1. Insert longitude and latitude of desired point in decimal °, save values in variables:
    point_Lon, point_Lat = transf(-1.23456789, 9.87654321)

These values must enter as numpy.ndarrays:

point_Lon = numpy.array([point_Lon])
point_Lat = numpy.array([point_Lat])

Now we can enter all this and get our results:

polarNeighs = wradlib.verify.PolarNeighbours(r=numpy.arange(1,nbins+1)*range_resolution,
                                             az=azi,
                                             sitecoords=[radar_Lon, radar_Lat],
                                             proj=projec,
                                             x=point_Lon, y=point_Lat,
                                             nnear=9)

From now on the tutorial (at step 6) can be continued normally.

Steps 1 and 2 are the whole point of my question, that escalated quickly!