BUG: Corrections converting gate to geographic location #353

carleyjmartin · 2023-08-16T22:19:59Z

Scope

This PR changes the calculation for geographic coordinates of the location of each range gate to include the 'offset' value of boresight from hdw files, and removes the abs() function from bmsep calculation to allow for the correct beam numbers.

As the abs() was originally a bug fix for a colouring-in issue with cartopy plots, this is now taken care of inside the fov method rather than in the position calculation using an if statement.

issue: #351

Approval

Number of approvals: 2 with a quick code review if possible

Test

matplotlib version: 3.7.1
Note testers: please indicate what version of matplotlib you are using

You can plot all the FOV in geographic coordinates, this should colour in the FOv in grey, anything outside should not be shaded:

import pydarn
from datetime import datetime
import matplotlib.pyplot as plt 

# North
stids = [209,208, 211,210,33,207,206,66,205,204,1,10,
         40,41,64,50,3,16,7,90,9,6,65,5,8,32]
plt.figure(figsize=(8, 8))
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5),
                        fov_color= 'grey',radar_location=True,
                        radar_label=True, coastline=True,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)

for stid in stids[1:]:
    _,_,ax,ccrs=pydarn.Fan.plot_fov(stid, datetime(2021, 2, 5, 0, 5), ax=ax,
                            ccrs=ccrs, fov_color= 'grey',
                            radar_location=True, radar_label=True,
                            coastline=True,
                            coords=pydarn.Coords.GEOGRAPHIC,
                            projs=pydarn.Projs.GEO)
plt.show()

# South
stids = [24,96,97,21,4,15,20,11,22,13,12,14,18,19]
plt.figure(figsize=(8, 8))
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5),
                        fov_color= 'grey',radar_location=True,
                        radar_label=True, coastline=True,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)

for stid in stids[1:]:
    _,_,ax,ccrs=pydarn.Fan.plot_fov(stid, datetime(2021, 2, 5, 0, 5), ax=ax,
                            ccrs=ccrs, fov_color= 'grey',
                            radar_location=True, radar_label=True,
                            coords=pydarn.Coords.GEOGRAPHIC,
                            projs=pydarn.Projs.GEO)
plt.show()

To check on the direction of beams, you can use the following code:

import pydarn
from datetime import datetime
import matplotlib.pyplot as plt 

plt.figure(figsize=(8, 8))
stids = [20]
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5),
                        fov_color= 'grey',radar_location=True,
                        radar_label=True, coastline=True,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5), ax=ax,
                        fov_color= 'red', beam=0, ccrs=ccrs,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5), ax=ax,
                        fov_color= 'blue', beam=15, ccrs=ccrs,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
plt.show()


plt.figure(figsize=(8, 8))
stids = [24]
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5),
                        fov_color= 'grey',radar_location=True,
                        radar_label=True, coastline=True,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5), ax=ax,
                        fov_color= 'red', beam=0, ccrs=ccrs,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
_,_,ax,ccrs=pydarn.Fan.plot_fov(stids[0], datetime(2021, 2, 5, 0, 5), ax=ax,
                        fov_color= 'blue', beam=15, ccrs=ccrs,
                        coords=pydarn.Coords.GEOGRAPHIC,
                        projs=pydarn.Projs.GEO)
plt.show()

This code checks the first and 15th beam (red and blue) for a -ve and +ve bmsep value:

These all require Cartopy to be installed - you can remove the projs, coords and coastline keywords to test in magnetic coordinates without cartopy.

Please check out the associated issue for further discussion!

egthomas · 2023-08-17T19:03:20Z

pydarn/utils/coordinates.py

-    beam_sep = np.radians(abs(SuperDARNRadars.
-                          radars[stid].hardware_info.beam_separation))
+    beam_sep = np.radians(SuperDARNRadars.
+                          radars[stid].hardware_info.beam_separation)
    rxrise = SuperDARNRadars.radars[stid].hardware_info.rx_rise_time

    # TODO: fix after slant range change


This if statement doesn't appear correct (which may explain the "TODO" note) - the beam_edge value should be zero if center is True (ie the current logic should be swapped).

Thanks Evan, I'll work through this and figure out what the proper logic we should be doing, there was a very long discussion about this two years ago, and I believe that it centred around what actually the definition of frang was (was it to the inner edge, or center of the first range gate) and the conclusion was essentially that no one actually had this definition set, different PIs thought it was different and it wasn't for us to define it at the time - as far as I remember, the outcome for us was that the logic for center was backwards from davitpy/rst but it made sense at the time. I'll link the main PR here, but there was a couple going on at the same time to figure this out: #134

@carleyjmartin yes, there still tends to be a lot of confusion around whether frang corresponds to the front edge or center of the first range gate (I think it's the latter). The convention for center vs edge of a beam in terms of azimuth should be more straightforward though: the center = True case should return the azimuth to the center of the beam. At least in the northern hemisphere, when center = False it should return the azimuth to the "left" / western" beam edge. Currently, it appears the logic in the python code is reversed from what it should be.

Ah, okay, yep. We have an issue here where we're using the variable name center in two similar but related ways. In range_separation the center value is in the direction from the radar outwards (which I assumed was the same thing in here!) and in coordinates is across the fov. I get it! I believe we may have flipped this one assuming the same definition of center being the same.

Right, so this is kinda confusing but, I think if we change the logic in the if statement around, we will then want to give center=False to the method from our plotting methods, as the plotting itself requires the values of the range gate corners.

So I think we also need to change the default to False. (EDIT: changing this in means we wont have to go and add the kw to the plotting methods that use this)

This makes more logical sense though, as if the user wants the central value they say center=True, otherwise we use False in general for plotting, and center retains the definition of giving the value at the center of the range gate.

I think we originally had it so that center was defined as True 'if the value we are given is the central value' which would cause the logic all to swap around and we had to then adjust the values to give us the edges to use for plotting. (if that makes sense, thinking of it from the 'what does the plot need' view rather than what logically does center mean)

EDIT: I've updated with this logic in the newest commit, I'll give it a proper test in the coming weeks. I just had a few hours free whilst it was raining at the site this morning!

egthomas · 2023-08-17T19:04:53Z

pydarn/utils/coordinates.py

@@ -172,7 +174,7 @@ def gate2geographic_location(stid: int, beam: int, height: float = None,
        beam_edge = 0

    # psi [rad] in the angle from the boresight
-    psi = beam_sep * (beam - offset) + beam_edge
+    psi = beam_sep * (beam - offset) - beam_edge + bmoff


Why was the sign changed from + beam_edge to - beam_edge ?

The sign change was suggested by Colin in this comment:
#351 (comment)

@carleyjmartin oh ok thanks, now I see that suggested change. This has the practical effect of shifting all of the beams eastward by one, eg printing out the psi values for the CVW radar for each beam:

0 -37.260000 -38.880000 -35.640000 1 -34.020000 -35.640000 -32.400000 2 -30.780000 -32.400000 -29.160000 3 -27.540000 -29.160000 -25.920000 4 -24.300000 -25.920000 -22.680000 5 -21.060000 -22.680000 -19.440000 6 -17.820000 -19.440000 -16.200000 7 -14.580000 -16.200000 -12.960000 8 -11.340000 -12.960000 -9.7200000 9 -8.1000000 -9.7200000 -6.4800000 10 -4.8600000 -6.4800000 -3.2400000 11 -1.6200000 -3.2400000 0.0000000 12 1.6200000 0.0000000 3.2400000 13 4.8600000 3.2400000 6.4800000 14 8.1000000 6.4800000 9.7200000 15 11.340000 9.7200000 12.960000 16 14.580000 12.960000 16.200000 17 17.820000 16.200000 19.440000 18 21.060000 19.440000 22.680000 19 24.300000 22.680000 25.920000 20 27.540000 25.920000 29.160000 21 30.780000 29.160000 32.400000 22 34.020000 32.400000 35.640000 23 37.260000 35.640000 38.880000

where the first column is beam number, the second column is for the beam center, the third column is using +beam_edge, and the fourth column is -beam_edge. You can see that the beam edges in the last column no longer correspond to the "left" beam edge, which is the standard convention. Maybe this is hemisphere specific?

Oooo okay, It could possibly be a fix for the -ve beam_sep valued radars, I think they're all in the southern hemisphere.

I'll figure it out!

@carleyjmartin I just did some tests in IDL, and I think it should still be + beam_edge for both Northern and Southern Hemisphere radars (as long as the beam_sep value retains the correct sign, which this pull request now ensures).

Thanks, Evan, I've switched it back around in the latest commit.

…lines

egthomas · 2023-08-25T13:40:39Z

@carleyjmartin thanks for making those changes - hopefully @Prof-CW agrees / can confirm!

carleyjmartin · 2023-09-05T19:14:15Z

Looks like @Prof-CW agrees with the changes over on the issue #351 , anyone feel free to merge if you like, or let me know if it needs any more changes! :)

carleyjmartin added 2 commits August 16, 2023 15:40

use bmoff and direction of bmsep

2be7d34

fix winding order in fov plots

ce51370

egthomas reviewed Aug 17, 2023

View reviewed changes

carleyjmartin and others added 2 commits August 21, 2023 10:20

rework for definition of center variable change, psi -ve to +ve, mod …

3e7795b

…lines

Merge branch 'develop' into ehn/psi-correction

9adbb66

carleyjmartin mentioned this pull request Sep 14, 2023

BUG: Update the position algorithm used for radar gate positioning #257

Closed

3 tasks

carleyjmartin mentioned this pull request Nov 13, 2023

RELEASE: pyDARN Version 4.0! #358

Merged

27 tasks

carleyjmartin merged commit 2a92b41 into develop Nov 23, 2023

carleyjmartin deleted the ehn/psi-correction branch November 23, 2023 16:32

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: Corrections converting gate to geographic location #353

BUG: Corrections converting gate to geographic location #353

carleyjmartin commented Aug 16, 2023

egthomas Aug 17, 2023

carleyjmartin Aug 21, 2023

egthomas Aug 21, 2023

carleyjmartin Aug 21, 2023

carleyjmartin Aug 21, 2023 •

edited

Loading

egthomas Aug 17, 2023

carleyjmartin Aug 21, 2023

egthomas Aug 21, 2023

carleyjmartin Aug 21, 2023

egthomas Aug 21, 2023

carleyjmartin Aug 21, 2023

egthomas commented Aug 25, 2023

carleyjmartin commented Sep 5, 2023

BUG: Corrections converting gate to geographic location #353

BUG: Corrections converting gate to geographic location #353

Conversation

carleyjmartin commented Aug 16, 2023

Scope

Approval

Test

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

carleyjmartin Aug 21, 2023 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

egthomas commented Aug 25, 2023

carleyjmartin commented Sep 5, 2023

carleyjmartin Aug 21, 2023 •

edited

Loading