Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

to_linecollection issues with projection and global grid wrapping w/ mpl and cartopy #984

Closed
zarzycki opened this issue Oct 1, 2024 · 4 comments · Fixed by #922
Closed
Labels
bug Something isn't working

Comments

@zarzycki
Copy link

zarzycki commented Oct 1, 2024

Version

v2023.08.4

How did you install UXarray?

Conda

What happened?

I am probably doing something unintelligent, but say I want to print the polylines assocated with a uxarray grid. In theory, this should work with mpl + cartopy...

import uxarray as ux
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

base_path = "./"
grid_path = base_path + "outCSne30.ug"
data_path = base_path + "outCSne30_vortex.nc"
uxds = ux.open_dataset(grid_path, data_path)

projection = ccrs.PlateCarree()  # also tried ccrs.Orthographic(central_longitude=0, central_latitude=90)
lc = uxds.uxgrid.to_linecollection(projection=projection,override=True)

#segments = lc.get_segments()
#print(segments)

fig, ax = plt.subplots(1, 1, figsize=(10, 10), constrained_layout=True,
                       subplot_kw={'projection': ccrs.Orthographic(central_longitude=0, central_latitude=90)})

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_collection(lc)
ax.set_global()
plt.title("Grid Orthographic Projection")
plt.show()

It gives me the projection but a blank graph:

Screenshot 2024-10-01 at 9 28 36 AM

I can sneak around this by creating the LineCollection with uxarray and then unpacking it and redefining it with a PC projection. The below now gives me grid lines, but it doesn't appear to wrap accurately along the date line (-180/180).

import uxarray as ux
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.collections import LineCollection

base_path = "./"
grid_path = base_path + "outCSne30.ug"
data_path = base_path + "outCSne30_vortex.nc"
uxds = ux.open_dataset(grid_path, data_path)

lc_original = uxds.uxgrid.to_linecollection()
segments = lc_original.get_segments()
projection = ccrs.PlateCarree()
lc = LineCollection(segments, linewidths=0.5, colors='black', transform=projection)

fig, ax = plt.subplots(1, 1, figsize=(10, 10), constrained_layout=True,
                       subplot_kw={'projection': ccrs.Orthographic(central_longitude=0, central_latitude=90)})

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_collection(lc)
ax.set_global()
plt.title("Grid Orthographic Projection")
plt.show()
Screenshot 2024-10-01 at 9 33 37 AM

So TL;DR two problems.

  1. Trouble getting to_linecollection() to be transformed correctly
  2. Lines not correctly wrapping around edges of mesh

If it's just something silly on my end (50/50 chance!), maybe we could add some additional lines to the mpl.ipynb notebook to help.

This may be tied to the periodic_elements option in to_polycollection() although the artifacts here aren't as pervasive.

What did you expect to happen?

Expected to use to_linecollection() to create a grid skeleton I could then draw with different projections using mpl + cartopy.

Can you provide a MCVE to repoduce the bug?

import uxarray as ux
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.collections import LineCollection

base_path = "./"
grid_path = base_path + "outCSne30.ug"
data_path = base_path + "outCSne30_vortex.nc"
uxds = ux.open_dataset(grid_path, data_path)

lc_original = uxds.uxgrid.to_linecollection()
segments = lc_original.get_segments()
projection = ccrs.PlateCarree()
lc = LineCollection(segments, linewidths=0.5, colors='black', transform=projection)

fig, ax = plt.subplots(1, 1, figsize=(10, 10), constrained_layout=True,
                       subplot_kw={'projection': ccrs.Orthographic(central_longitude=0, central_latitude=90)})

ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_collection(lc)
ax.set_global()
plt.title("Grid Orthographic Projection")
plt.show()
@zarzycki zarzycki added the bug Something isn't working label Oct 1, 2024
@philipc2
Copy link
Member

philipc2 commented Oct 1, 2024

Hi @zarzycki

The below now gives me grid lines, but it doesn't appear to wrap accurately along the date line (-180/180).

By default, the values along the date line are being excluded, since they require additional processing to split when we convert our data into geometries. This is indicated by the periodic_elements parameter.

Here's a reference.

# split periodic elements (i.e. those along the dateline)
lc_original = uxds.uxgrid.to_linecollection(periodic_elements='split')

This works with your work around.

image

There appears to be a bug in the projection code. I'll get this fixed as part of #922

@zarzycki
Copy link
Author

zarzycki commented Oct 1, 2024

Thanks @philipc2!

The user API doesn't mention periodic_elements: https://uxarray.readthedocs.io/en/v2024.02.0/user_api/generated/uxarray.Grid.to_linecollection.html

Although I thought I tried it and got an error. Maybe I messed something up...

FWIW, while that plot looks better, there does appear to be a weird collection of line segments right at the pole. Any idea if that's the mesh or the plotting? I can look later when I have some free time, too.

@philipc2
Copy link
Member

philipc2 commented Oct 1, 2024

https://uxarray.readthedocs.io/en/v2024.02.0/user_api/generated/uxarray.Grid.to_linecollection.html

That appears to be an older release. here's our latest documentation:

https://uxarray.readthedocs.io/en/v2024.08.2/user_api/generated/uxarray.Grid.to_linecollection.html

FWIW, while that plot looks better, there does appear to be a weird collection of line segments right at the pole. Any idea if that's the mesh or the plotting? I can look later when I have some free time, too.

That's most likely a slight artifact of the split polygons. I can look into it further.

@zarzycki
Copy link
Author

zarzycki commented Oct 1, 2024

Gah, I blame Google but I should have checked that.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
Status: ✅ Done
Development

Successfully merging a pull request may close this issue.

2 participants