SuperDARN / pydarn

Python library for visualizing SuperDARN Data
GNU Lesser General Public License v3.0
31 stars 11 forks source link

psi calculation for beam,range coords #351

Closed Prof-CW closed 9 months ago

Prof-CW commented 1 year ago

Discussion

Calculating psi from beam_sep in function gate2geographic_location

Description

In /pydarn/utils/coordinates.py and function gate2geographic_location, psi is computed as: psi = beam_sep * (beam - offset) + beam_edge

Just above that we have the following code snippet; !# Some beam seperations are negative which changes how the coordinates wrap !# we absolute to make it easier of fov-color for cartopy plotting beam_sep = np.radians(abs(SuperDARNRadars.radars[stid].hardware_info.beam_separation))

so that beam_sep is always +ve.

For example, with this code, it looks like psi for beam 0 puts it on the east rather than west side of the boresite.

Am interested to know if eastward scanning radars (start from the west) need to preserve the -ve beam_sep compared with westward scanners (that start in the east).

Category

Documentation

This is a follow on from #257

egthomas commented 1 year ago

This is tangentially related, but there appears to be a missing bmoff term in the psi calculation (which accounts for electronic shifts to the radar boresite), eg

    /* Calculate deviation from boresite [deg] */
    psi=pos->bmsep*(bcrd-offset)+bm_edge+pos->bmoff;

Note the bmoff parameter was only added to the hardware files ~2 years ago. I think some examples of radars with time-dependent, non-zero bmoff values were BKS, WAL, and one of the Syowa radars.

carleyjmartin commented 1 year ago

All of these questions are really making me dig into my extremely fallible memory haha! We did have some issues a while back with winding order (colors everything outside the fov instead of inside) for the beams and gates, and I believe the 'abs' may have been the fix for that in the beam_sep line (I went back in the blame but this section has been bounced around between modules so I can't quite figure it's original commit right now).

There is a new feature recently put into the develop branch that plots specific beams in a fan plot so we can use that to make sure its plotting the beams from the correct direction #298 . I'm happy to do this, I just won't be able to do it for a few weeks as I've got some other commitments.

With regards to the equation with 'bmoff' in @egthomas, I'm sorry I very rarely use C, so I'm not really sure what that line is programatically doing? Otherwise, I'm also happy to add the bmoff in, especially if it'll be effecting some of the radar beam positions!

egthomas commented 1 year ago

With regards to the equation with 'bmoff' in @egthomas, I'm sorry I very rarely use C, so I'm not really sure what that line is programatically doing? Otherwise, I'm also happy to add the bmoff in, especially if it'll be effecting some of the radar beam positions!

@carleyjmartin no problem, I should have provided a little more detail in my comment. Currently, it appears the psi calculation in gate2geographic_location is

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

whereas in the C and IDL code in the RST, we have added support for one additional term (bmoff) which accounts for an additional shift in the radar boresight direction which is tracked in the radar hardware files, ie

    /* Calculate deviation from boresite [deg] */
    psi = pos->bmsep * (bcrd - offset) + bm_edge + pos->bmoff;

Note this psi calculation is nearly identical to the pydarn code, except for some different variable names and that missing bmoff term (also note I've added some spaces for readability).

Here is an excerpt from the WAL hardware file showing the electronic offset (ie bmoff) of -9.8 degrees which must be applied to the radar beams / field of view between Nov 2005 and Mar 2006 (OFF SET column):

#                            GEOGRAPHIC             BORE-   OFF  BEAM V  |--------- interferometer ---------| |--- RX ---| RNG MX
#      YYYYMMDD HH:MM:SS    LAT      LON      ALT   SIGHT   SET  SEP  SG PS  TDIFF   TDIFF   X      Y     Z    RXRIS  ATTN     BM
 32  1 20050610 00:00:00   37.930  -75.470    50.0   35.9  0.00  1.62  1  1  0.000   0.000   0.0  100.0   0.0    0.0   0 0 110 16
 32  1 20051104 18:00:00   37.930  -75.470    50.0   35.9 -9.80  1.62  1  1  0.000   0.000   0.0  100.0   0.0    0.0   0 0 110 16
 32  1 20060101 12:00:00   37.930  -75.470    50.0   35.9 -9.80  3.24  1  1  0.000   0.000   0.0  100.0   0.0    0.0   0 0 110 16
 32  1 20060313 19:38:00   37.930  -75.470    50.0   35.9  0.00  3.24  1  1  0.000   0.000   0.0  100.0   0.0    0.0   0 0 110 24
carleyjmartin commented 1 year ago

Thanks Evan! I was just wondering what the -> meant in C but I think it's the equivalent of pos['bmoff'] in python, where pos is an object containing the hdw data? So, it's just added to the end 🙈

Prof-CW commented 1 year ago

Of course, the hdw file lists values in degress, but I think the python code does the calculation in radian so bmoff would need to be in rad.

For the GEO locations of the data, I checked where beam0 data were being mapped in GEO using (Basemap) for south hemisphere using

f_geo_lats, f_geo_lons = pydarn.Coords.GEOGRAPHIC(hdw_dat.stid,date=dt.datetime(yy,mnth,dd), beams=15,gates=(0,75), range_estimation=pydarn.RangeEstimation.SLANT_RANGE) geo_lats, geo_lons = pydarn.Coords.GEOGRAPHIC(hdw_dat.stid,date=dt.datetime(yy,mnth,dd), beams=0,gates=(0,75), range_estimation=pydarn.RangeEstimation.SLANT_RANGE)

f_x, f_y = m(f_geo_lons,f_geo_lats) m.scatter(f_x,f_y,s=0.5,color='grey') x, y = m(geo_lons,geo_lats) m.scatter(x,y,s=0.5,color='black')

to give beam 0 mapping (black) and other beams (grey) as: sth_fov

I know that beam 0 for both BPK and TIG are on the wrong side of the boresite here. Beam 0 for BPK and TIG should map to the west side of the boresite. I am not sure about the others. For all 5 radars mapped here, the hdw files have hdw_dat.beam_separation -ve, which is being removed with the abs on beam_sep in coordinates.py.

If I change two lines in gate2geographic_location (in coordinates.py) as follows: beam_sep = np.radians(SuperDARNRadars.radars[stid].hardware_info.beam_separation)

i.e. remove the abs and change the sign on beam_edge to

psi = beam_sep * (beam - offset) - beam_edge

then I get the beam 0 mapping on the correct side (for BPK and TIG) as

sth_fov_clw

The beam scan direction for BPK is east to west, but the beam sequence for this is high beam number to low beam number with beam0 the western-most beam as shown in the modified map.

carleyjmartin commented 1 year ago

Hi both, I've amended the lines with your suggestions so that they now read:

beam_sep = np.radians(SuperDARNRadars.
                          radars[stid].hardware_info.beam_separation)

psi = beam_sep * (beam - offset) - beam_edge + bmoff

If I plot each of the southern hemisphere radars on a fov plot, highlighting beam 0 in red on each, I think that it appears to be correct for all the -ve bmsep radars? The plot may look a little weird compared to yours above, this is plotted from above the north pole through the earth in AACGMv2.

Screenshot 2023-08-16 at 3 37 14 PM

I'm not sure if this is correct or not but we have been reading in the offset value in the hdw files as a secondary boresight called 'electronic' (the regular boresight is called 'physical'). Are these descriptions correct?

I'm going to have to fiddle with some of the plotting though to stop the winding order colouring outside the FOV instead of inside, so it's not quite ready for a pull request yet. The code can be found on the branch called ehn/psi-correction which I will push to in a couple of minutes. - was an easier fix than I expected, there should be a link directly below to the PR :) Thanks!

Prof-CW commented 1 year ago

@carleyjmartin Well, the plot looks weird to me but I understand why you did it that way, with east and west switched and looking through the Earth from the north. I have met a few who have trouble with mapping in the south hemisphere! Beam 0 for BPK, UNW and TIG look right, being on the west side of the boresite and scanning eastward. I have not checked the other radars. Might be a question/confirmation from the relevant PIs for those.

Prof-CW commented 1 year ago

Looks good for removing the abs(). The only item to check is whether the psi calculation should be; psi = beam_sep (beam - offset) - beam_edge + bmoff or psi = beam_sep (beam - offset) + beam_edge + bmoff These differ by one beam width in the FOV span. It looks like Rob Barnes had the + beam_edge version according to @egthomas so we had probably best do that, rather than the -ve beam_edge that @carleyjmartin has above.