SuperDARN / pydarn

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

Faster coastline coordinate conversion #380

Closed billetd closed 5 months ago

billetd commented 5 months ago

Scope

This PR changes geo_coastline_to_mag() so it works on geometry arrays, rather than iterating over each glat/glon pair individually, when converting to AACGM. The result is a significant speedup.

Approval

Number of approvals: 2

Test

Taken from https://pydarn.readthedocs.io/en/main/user/fan/, but can be used in any projection where coastline=True and Cartopy is installed.


import pydarn

file = "path/to/fitacf/file"
SDarn_read = pydarn.SuperDARNRead(file)
fitacf_data = SDarn_read.read_fitacf()

pydarn.Fan.plot_fan(fitacf_data,
                    scan_index=1, lowlat=60, zmin=-1000, zmax=1000,
                    boundary=True, radar_label=True,
                    groundscatter=True, ball_and_stick=True, len_factor=300,
                    coastline=True, parameter="v")
plt.show()
carleyjmartin commented 5 months ago

So this works really well, but I appear to have lost mainland Russia? I'm guessing maybe an indexing issue when going through the geoms?

Screenshot 2024-03-26 at 12 06 18 PM

EDIT: Russia is a multi-part geometry that is breaking int he try and not plotting. I ahve code somewhere that I'll find that amends this.

carleyjmartin commented 5 months ago

This is how I dealt with multi-part geometries before, probably have to do something similar in the projections file:

import cartopy.feature as cfeature
import datetime as dt
from pydarn.plotting.projections import convert_geo_coastline_to_mag
from shapely.geometry import mapping

# Set a time, required for the magnetic conversion
time = dt.datetime(2024, 1, 18, 0 ,0)

# Read in the geometry object of the coastlines
cc = cfeature.NaturalEarthFeature('physical', 'coastline', '50m',
                                      color='k', zorder=2.0)

# Iterate over the geometry objects (with coordinates in geographic) and send them to get converted
for geom in cc.geometries():
    if geom.__class__.__name__ == 'MultiLineString':
        for g in geom.geoms:
            coastline_mlt = convert_geo_coastline_to_mag(g, time)
            coastline_mlon = convert_geo_coastline_to_mag(g, time, mag_lon=True)
            mlt, mlat = coastline_mlt.coords.xy[0], coastline_mlt.coords.xy[1]
            mlon, mlat = coastline_mlon.coords.xy[0], coastline_mlon.coords.xy[1]
    else:
        coastline_mlt = convert_geo_coastline_to_mag(geom, time)
        coastline_mlon = convert_geo_coastline_to_mag(geom, time, mag_lon=True)
        mlt, mlat = coastline_mlt.coords.xy[0], coastline_mlt.coords.xy[1]
        mlon, mlat = coastline_mlon.coords.xy[0], coastline_mlon.coords.xy[1]
carleyjmartin commented 5 months ago

Issue was with how we shifted longitudes, if the first longitude was a NaN it used that to shift all in the geom resulting in all NaNs.

Also added the option to specify the cartopy_scale for '10m', '50m' and '110m', 110 is default, 50 and 10 will slow down your plotting considerably so be aware. In doing this I also fixed the issue with multipart geometires being plotted, so it can convert 50m and 10m scales to geomag.

Please test for maps and fan at least in geomag and geographic coordinates. You can use the key word 'cartopy_scale='50m'' to test the scales on anything that you're plotting a coastlines with.

Doreban commented 5 months ago

Using matplotlib version 3.8.3

I tested making maps and fan plots using the new branch for cartopy scales of '10m', '50m', and '110m'. Plots came out as expected, all the coastlines expected were there, and the changing of the cartopy scale did increase the resolution of the coastlines on the plots.

Doreban commented 5 months ago

I also tested both geomag and geographic coordinate systems for the fan plots and it worked.