BUG: Update the position algorithm used for radar gate positioning

[carleyjmartin]

BUG

This is more of a question or a sanity check for me. Currently, we have a default 180 km set as the distance to the first range gate frang value if no value is read in or given. However, when we call lat,lon=pydarn.Coords.GEOGRAPHIC(5) to get the location of range gates from the Saskatoon radar, the actual lat/lon returned is around 46 km away from the radar and not 180 km, I may expect some slight differences as we're doing a number of geodetic transformations on this, but I'm not sure that this large a difference is true?

Priority

• [ ] Minor - typo, incorrect naming, fails on specific uncommon situation
• [x] Major - misrepresentation of data, failure in a common situation
• [ ] Urgent - failure in usage, large misrepresentation

Information

python version: 3.9 OS: MacOS 10.15.5 Catalina matplotlib version: 3.5.1

Example of the bug

FOV = pydarn.Coords.GEOGRAPHIC(5)
lat = FOV[0]
lon = FOV[1]

Attempts

What have you tried already to fix it? I have made no attempt to fix, other than search for the problem. I believe that it is reading in the slant ranges from range_estimation correctly, but somewhere in geocentric_coordinates we're loosing a few kms.

Data Location

Does not require data.

Potential Bug Location

Possibly in geodetic2geocentric in geo.py

[carleyjmartin]

It appears that the algorithms used in geo.py to calculate the position of range gates may be out of date/better may be available. Requires more investigation when someone has time.

[ecbland]

This is the code used in RST for determining the geographic and geomagnetic coordinates of the FOV cells: https://github.com/SuperDARN/rst/blob/develop/codebase/superdarn/src.lib/tk/rpos.1.7/src/cnvtcoord.c There might be some differences with the rbposlib IDL code that was used for pyDARN

[Prof-CW]

I came across this 'bug' when recently looking at the ground scatter coord mapping. I keep getting an error (raised by a -ve sqrt) when I try to get the km to a range gate assuming ground scatter. I hunted this down to a function (gate2groundscatter) that is defined in the range_estimations.py code (modified aug 2022). The comment at the top of 'gate2groundscatter' says the eqn used is based on the Bristow et al. paper (1994) pg 325. Well I looked at the paper, looked at that eqn, drew a lot of range-distance pictures, and found that the eqn used is quite the approximation. The variable, ground_scatter_mapped_ranges, has the sqrt before the arcsin and one can easily see why it is going -ve and popping up an error. I do have an alternative solution (eqn) but it requires elevation angle is known. Before proceeding, thought I would check to see if any progress has been made on this and if not, perhaps can work it together. Let me know ColinW (Newcastle Uni, Australia)

[egthomas]

@Prof-CW I agree about the Bristow et al equation - I raised similar concerns about its validity a few times in the past (eg, https://github.com/SuperDARN/pydarn/pull/183 and https://github.com/vtsuperdarn/davitpy/issues/353). For a more generalized set of geolocation routines, you can see the rpos_v2 functions we recently added to the RST:

https://github.com/SuperDARN/rst/blob/main/codebase/superdarn/src.lib/tk/rpos.1.7/src/rpos_v2.c

as well as a simple program to test them:

https://superdarn.github.io/rst/superdarn/src.bin/tk/testing/rpos_v2_test/index.html

[w2naf]

Hi @egthomas,

You wrote in #183 that the Bristow equation is solving for the ground range to the ionospheric reflection point (D) and not the slant range to the reflection point (R/2). I think we technically don't want either of those for the type of calculations we are doing with TIDs. I think the thing we really want is range at reflection height from the radar to the reflection point. Or in other words, the ground range to the reflection point projected up to ionospheric heights.

[egthomas]

@w2naf thanks for the reply. You're right that the first part of my comment has since been addressed (ie the slant range vs ground range calculation). The latter (and I believe still relevant) concern was

I'm also not convinced the equation in that paper is even correct for solving the ground range to the ionospheric reflection point of a ground scatter echo, but that's another matter.

The attached plot shows the ground range to the ionospheric reflection point of 1-hop ground scatter calculated using the Bristow equation (black), a flat Earth assumption (red), and a spherical Earth assumption (blue); the bottom panel shows the differences between each method. Granted the differences are not large relative to the typical range gate separation (45 km), but they do exist. [edit: note this is assuming a 300 km reflection height at all slant ranges]

Maybe I'm just misunderstanding what you mean by the ground range to the reflection point projected up to ionospheric heights? How does that differ from the ground range of the ionospheric reflection point?

[carleyjmartin]

On the non-scientific coding side (and a definite non-expert on geolocation):

More than happy to have another alternative range estimation for users, fairly easy to implement all this, @Prof-CW is there an article which we can link for a new equation?

This issue was mainly about the 'rbpos' code that was used to write the geolocation algs was not the most up to date, I believe we need to update to the 'rpos_v2' that @egthomas linked above (thank you!)

If we come to a consensus on what we need, I'm happy to implement, I'll just need to be told what to do! Very happy we have some traction on this issue as it's been bugging (pun intended) me for a long time!

If you'd like a forum to discuss this, our next pyDARN meeting I believe is next Friday (21st July) at 18 UTC but happy to move time to make attendance easier if that's okay with @bharatreddy

[Prof-CW]

Thanks for some of the history on this one. I have also noticed that the lat,lons (both GEO and AACGM) output from pydarn are different to those output from the VT website under "Data Diagnostics" -> "Print Fit Records". This started me looking at the calculation in pydarn. Interesting to also learn that rpos has been added to RST. One would think that the lat,lon outputs from RST and pydarn should be the same for all radars. I am going back to the maths and diagrams and working forward. That way we should get both RST and pydarn correct (and the same) with the assumptions clear (flat Earth, curved Earth etc.).

It is relatively easy to derive the Matlab slant range equations by applying the sine and cosine rules to the ray trace diagram. See https://www.mathworks.com/help/radar/ref/slant2range.html under "More About". These give equations that require elevation angle. Has anyone derived the equations used in rbpos_v2 ? If not, am happy to have a go at it. C

[w2naf]

Hi @egthomas,

How is the spherical-Earth ground range that you show calculated? I would have thought that would be the same as the Bristow equation.

When I say "ground range to the reflection point projected up to ionospheric heights," I mean length D in Figure 2b of https://doi.org/10.1002/2014JA019870 projected up to the height where the TID is occurring (say 300 km altitude).

Thanks for your help with this.

[Prof-CW]

Yes, I think that is right. I am using this diagram where h is virtual height and D is range along the ground (circular path). Angle {epsilon} is elevation angle. Assume we do not need to add an antenna height off the ground.

Seems to me the Bristow et al. eqn (eqn 1 in https://doi.org/10.1002/2014JA019870) is from: D = Re alpha (arc length along the ground) d = sqrt( r^2 - h^2) (pythagorus of right-angle triangle, my r is your R/2) and sin (alpha) = d/Re (sine definition from right angle triangle) The approximation comes from the relationship between D, d, and h. For D not close to d, h is less than it should be.

After we get D, does anyone have a justification for using 'great circle' distance rather than using the distance in a plane formed by the 2 points defining r and the Earth centre?

[Prof-CW]

I have attempted to collect slant and ground range calculation equations and justifications into one document. This might addess @carleyjmartin question about papers that might be linked to the relevant equations. I have also checked the equations in the RST C code that @egthomas linked and they appear to be correct subject to the assumptions I have listed in the pdf. rpos.pdf

@carleyjmartin let me know if you would like me to check the pydarn code for these calculations. I suggest we swap out the ground scatter Bristow approximation (Eqn 1 in pdf) in pydarn to be the same as the RST C code which is, Eqn 11 in the pdf.

[egthomas]

@Prof-CW you may also want to see our recent virtual height model paper where we described (most) of the equations in the RST C code I linked:

https://doi.org/10.1029/2022RS007429

[carleyjmartin]

Hi all,

I had a go at adding in a new range estimation with the new equation (11). I've kept the old version in just in case but now it's called RangeEstimation.GSMR_BRISTOW, with the newer equation being the RangeEstimation.GSMR.

I noticed that this new equation is much closer to the general slant range values than the old GSMR calculation is, whereas using the Bristow equation is much closer to values given from half of slant range. Is this expected, or am I missing a /2 somewhere?

The new equation is on branch ehn/new-gsmr-eqn, using the following code you can see the difference in the range time plots below:

import matplotlib.pyplot as plt 
import pydarn
file = "/Users/carley/Documents/data/20211102.0000.00.rkn.a.fitacf"
SDarn_read = pydarn.SuperDARNRead(file)
fitacf_data = SDarn_read.read_fitacf()
a = pydarn.RTP.plot_summary(fitacf_data, beam_num=2, range_estimation=pydarn.RangeEstimation.GSMR, watermark=False, title='New GSMR')
plt.show()
a = pydarn.RTP.plot_summary(fitacf_data, beam_num=2, range_estimation=pydarn.RangeEstimation.GSMR_BRISTOW, watermark=False, title='Old GSMR')
plt.show()
a = pydarn.RTP.plot_summary(fitacf_data, beam_num=2, range_estimation=pydarn.RangeEstimation.SLANT_RANGE, watermark=False, title='Slant Range')
plt.show()

And just so you don't have to go digging, this is the new code in utils/range_estimations.py. The old gate2groundscatter() is now gate2gs_bristow() and gate2groundscatter() contains the new equation now.

def gate2gs_bristow(reflection_height: float = 250, **kwargs):
    """
    Calculate the ground scatter mapped range (km) for each slanted range
    for SuperDARN data. This function is based on the Ground Scatter equation
    from Bristow paper at https://doi.org/10.1029/93JA01470 on page 325
    Parameters
    ----------
        reflection_height: float
            reflection height
            default:  250

    Returns
    -------
        ground_scatter_mapped_ranges : np.array
            returns an array of ground scatter mapped ranges for the radar
    """
    slant_ranges = gate2slant(**kwargs)
    ground_scatter_mapped_ranges =\
            Re*np.arcsin(np.sqrt((slant_ranges**2/4)-
                                 (reflection_height**2))/Re)

    return ground_scatter_mapped_ranges

def gate2groundscatter(reflection_height: float = 250, hop: float = 0.5,
                       **kwargs):
    """
    Calculate the ground scatter mapped range (km) for each slanted range
    for SuperDARN data. This function is based on the Ground Scatter equation
    discussed in the issue github.com/SuperDARN/pydarn/issues/257
    Parameters
    ----------
        reflection_height: float
            reflection height
            default:  250
        hop: float
            hop numberof returning data
            default: 0.5

    Returns
    -------
        ground_scatter_mapped_ranges : np.array
            returns an array of ground scatter mapped ranges for the radar
    """
    slant_ranges = gate2slant(**kwargs)

    num = - Re**2 - (Re + reflection_height)**2\
          + (slant_ranges * (0.5 / hop))**2
    den = 2 * Re * (Re + reflection_height)
    ground_scatter_mapped_ranges = (hop/0.5) * Re * np.arccos( - num / den )

    return ground_scatter_mapped_ranges

I've included the option to add in that 'hop' variable with a default of 0.5, but for now I don't really know how to properly implement that if it changes, unless the user sets it.

Let me know if you see anything suspicious!

I'm pretty sure we (I) need to go in and have a proper look at the geo.py module that contains the code to convert from these range estimations to lat and lon values in geo and mag coordinates, as this is where I'm seeing some issues with that first range gate distance in all the range estimations. But that can be a different pull request when I get around to it.

[Prof-CW]

@carleyjmartin Looks like you need a divide by 2 in the gatetogroundscatter function, i.e. num = - Re**2 - (Re + reflection\_height)**2 + (slant\_ranges/2 * (0.5 / hop))**2 The old (Bristow) and revised eqns should give similar ground ranges (D) for the closer slant ranges as per the graphs by @egthomas posted 13 July. The /2 will depend on whether your slant_ranges variable is r or 2r in Figure 1 of the pdf (and posted above). Your gate2gs_bristow function looks like it already has this slant_ranges/2 since you have: ground_scatter_mapped_ranges = Re*np.arcsin(np.sqrt((slant_ranges**2/4-(reflection_height**2))/Re) The slant_ranges**2/4 in there is (slant_ranges/2)**2 Also suggest we add a trap for -ve sqrt in the bristow function, which can happen for close slant ranges and ill informed virtual height inputs. We should also be clear in the comments, replacing reflection height with virtual height.

I will head to geo.py and take a look at the ground range to lat lon mapping.

[carleyjmartin]

Thanks @Prof-CW I have included the /2 and renamed the reflection height to virtual height in both methods. The new method using eqn 11 looks far more similar to the old method using eqn 1 now.

For the sqrt issue: We need to retain the shape of the array [#beams, #rangegates] to return to all the plotting methods and keep them happy. These include the inf values and then the plotting methods deal with the inf values there. To retain this shape, I have added a check in the Bristow method to see if there are inf values, and if there are then we show a little warning but the plot will still go ahead and plot, usually just by removing the inf values at plotting time.

There is currently a bug with warnings in summary plots that has a fix still to be tested in #328 so I have tested this warning in range time plots alone. Just in case anyone tests and sees that the warning is not coming up - it will once we merge the other PR I promise!

This is the newest version of the code:

def gate2gs_bristow(virtual_height: float = 250, **kwargs):
    """
    Calculate the ground scatter mapped range (km) for each slanted range
    for SuperDARN data. This function is based on the Ground Scatter equation
    from Bristow paper at https://doi.org/10.1029/93JA01470 on page 325
    Parameters
    ----------
        virtual_height: float
            virtual height
            default:  250

    Returns
    -------
        ground_scatter_mapped_ranges : np.array
            returns an array of ground scatter mapped ranges for the radar
    """
    slant_ranges = gate2slant(**kwargs)
    ground_scatter_mapped_ranges = Re * np.arcsin(np.sqrt((slant_ranges**2 / 4)
                                                  - (virtual_height**2)) / Re)
    # Check to see if there is an issue with the sqrt and
    # give user a warning if so, these values will be dealt with in
    # the individual plotting algs as we need to return the full array
    # of values for the complete beam*range gate array
    if any(np.isfinite(ground_scatter_mapped_ranges)):
        warnings.warn("Warning: Be aware that the range estimation"
                      " you have chosen has calculated some infinite"
                      " values. These values will not be plotted."
                      " You may use RangeEstimation.GSMR to avoid this"
                      " issue.")
    return ground_scatter_mapped_ranges

def gate2groundscatter(virtual_height: float = 250, hop: float = 0.5,
                       **kwargs):
    """
    Calculate the ground scatter mapped range (km) for each slanted range
    for SuperDARN data. This function is based on the Ground Scatter equation
    discussed in the issue github.com/SuperDARN/pydarn/issues/257
    Parameters
    ----------
        virtual_height: float
            virtual height
            default:  250
        hop: float
            hop number of returning data
            default: 0.5

    Returns
    -------
        ground_scatter_mapped_ranges : np.array
            returns an array of ground scatter mapped ranges for the radar
    """
    slant_ranges = gate2slant(**kwargs)

    num = - Re**2 - (Re + virtual_height)**2\
        + (slant_ranges/2 * (0.5 / hop))**2
    den = 2 * Re * (Re + virtual_height)
    ground_scatter_mapped_ranges = (hop/0.5) * Re * np.arccos(- num / den)

    return ground_scatter_mapped_ranges

[carleyjmartin]

Thanks for the discussion everyone, so far we have two PRs open to test relating to this issue: #356 and #353