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

Support GEOS Cube-Sphere Grids #802

Merged
merged 12 commits into from
Jun 25, 2024
Merged

Support GEOS Cube-Sphere Grids #802

merged 12 commits into from
Jun 25, 2024

Conversation

philipc2
Copy link
Member

@philipc2 philipc2 commented May 24, 2024

Closes #801

Overview

  • Add support for reading and encoding the GEOS Cube-Sphere format in the UGRID conventions

@philipc2
Copy link
Member Author

Plot of a ~200,000,000 element grid (CS5760)

image

@philipc2
Copy link
Member Author

Plot of a ~900 element grid (CS12)

image

@philipc2 philipc2 changed the title DRAFT: Support GEOS Cube-Sphere Grids Support GEOS Cube-Sphere Grids May 29, 2024
@philipc2 philipc2 marked this pull request as ready for review May 29, 2024 20:24
@philipc2
Copy link
Member Author

@mathomp4

Let me know if you have any feedback!

@mathomp4
Copy link

@mathomp4

Let me know if you have any feedback!

Well, I think I need to learn more Python but I decided to try on a stretched grid file I finally got. I first tried plotting the grid but apparently this:

<uxarray.Grid>
Original Grid Type: GEOS-CS
Grid Dimensions:
  * n_node: 440646
  * n_face: 437400
  * n_max_face_nodes: 4
  * n_nodes_per_face: (437400,)
Grid Coordinates (Spherical):
  * node_lon: (440646,)
  * node_lat: (440646,)
  * face_lon: (437400,)
  * face_lat: (437400,)
Grid Coordinates (Cartesian):
Grid Connectivity Variables:
  * face_node_connectivity: (437400, 4)
Grid Descriptor Variables:

was a bit too much for my laptop. I waited a few minutes and in the end it was just solid blue.

But following some of your examples, I did get out:

c270_plot

which is nice! I then tried a plot.polygons() and my computer hated me again. 😄

I think what this demonstrates is that I need to convince my boss to help make a much coarser stretch grid...or learn more UXarray

@philipc2
Copy link
Member Author

Glad you were able to get it working! Surprised that the Polygon plot took that long. Typically grids of that size should take no longer than a minute maximum for the first run.

@mathomp4
Copy link

Glad you were able to get it working! Surprised that the Polygon plot took that long. Typically grids of that size should take no longer than a minute maximum for the first run.

Hmm. Well, it was on an M3 MacBook Air...if that matters?

Note: it's also possible I called things incorrectly. I was just sort of trying to echo the cookbook. 😄

uxarray/io/utils.py Show resolved Hide resolved
@philipc2
Copy link
Member Author

Glad you were able to get it working! Surprised that the Polygon plot took that long. Typically grids of that size should take no longer than a minute maximum for the first run.

Hmm. Well, it was on an M3 MacBook Air...if that matters?

Note: it's also possible I called things incorrectly. I was just sort of trying to echo the cookbook. 😄

Yeah, I've run similar sized meshes on a 32GB M1 and had reasonable performance. Mind sharing the file so I can try to recreate it?

@mathomp4
Copy link

Yeah, I've run similar sized meshes on a 32GB M1 and had reasonable performance. Mind sharing the file so I can try to recreate it?

Yes...if I can figure out how to share it. Can you try:

https://gmao.gsfc.nasa.gov/gmaoftp/gmao_SIteam/test-fixapp-c270.geosgcm_prog_nat.20200415_0000z.nc4

I think that works. If you get it, let me know. It's about 500 MB because I didn't really clean it up

@philipc2
Copy link
Member Author

import uxarray as ux
import geoviews.feature as gf
import cartopy.crs as ccrs

features = gf.coastline(
    projection=ccrs.PlateCarree(), line_width=1, scale="50m"
) * gf.states(projection=ccrs.PlateCarree(), line_width=1, scale="50m")

file_path = r"C:\Users\chmie\Downloads\test-fixapp-c270.geosgcm_prog_nat.20200415_0000z.nc4"
%%time
uxds = ux.open_dataset(file_path, file_path)

Took about 3s

%%time
uxds['T'][0][-1].plot.rasterize(method='polygon', exclude_antimeridian=True) * features

Took about 3s as well.

image

uxds['RH'][0][0].subset.bounding_circle((-90, 25), 10).plot.rasterize(method='polygon', exclude_antimeridian=True, pixel_ratio=4.0) * features

image

@philipc2
Copy link
Member Author

I'd like to mention that the white spots are due to an upstream issue with Datashader, which will be fixed in the next release.

holoviz/datashader#1329

@mathomp4
Copy link

Hmm. When I first tried your first block of code, I got:

ValueError: No plotting extension is currently loaded. Ensure you load an plotting extension with hv.extension or import it explicitly from holoviews.plotting before applying any options.

So, I added:

import holoviews as hv
hv.extension("bokeh")

And then the plots took about 2s. So...yeah, I must have really screwed things up in my testing to be ~100x slower. 😄

@mathomp4
Copy link

And I was wondering about the weird white spots, but I figured they were due to just...something.

Now, my laptop was not happy with:

uxgrid = ux.open_grid(file_path,file_path)

%%time
uxgrid.plot()

Actually, my browser sort of crashed before it could render...

@philipc2
Copy link
Member Author

uxgrid = ux.open_grid(file_path,file_path)

%%time
uxgrid.plot()

And I was wondering about the weird white spots, but I figured they were due to just...something.

Now, my laptop was not happy with:

uxgrid = ux.open_grid(file_path,file_path)

%%time
uxgrid.plot()

Actually, my browser sort of crashed before it could render...

The issue here is that you're using open_grid instead of open_dataset. The default plotting for a Grid object is an edge plot of the mesh, which is indeed very slow for any mesh larger than a few thousand elements.

@mathomp4
Copy link

Ohhhh. So:

uxds = ux.open_dataset()
uxds.uxgrid.plot()

would be faster? I'll try that tonight on my laptop!

(I'm learning!)

@philipc2
Copy link
Member Author

Ohhhh. So:

uxds = ux.open_dataset()
uxds.uxgrid.plot()

would be faster? I'll try that tonight on my laptop!

(I'm learning!)

Plotting a data variable would be faster, since that performs a raster operation as opposed to a vector plot of the geometries.

Our grid (i.e. only mesh) plotting is all vector based, but our data variables can be plotted as rasters.

# defaults to raster plot
uxds['T'].plot()

# can specifically call a raster
uxds['T'].plot.rasterize(method='polygon')

Copy link
Member

@aaronzedwick aaronzedwick left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice work on this!

@mathomp4
Copy link

Nice work on this!

Doubly so from me and the GEOS community. Nice to have another package that can read our files without need to do all sort of gymnastics!

@philipc2 philipc2 merged commit 13462bd into main Jun 25, 2024
15 checks passed
@mathomp4
Copy link

Woo! Thanks. I eagerly await the next release! 😄

@mathomp4
Copy link

I'd like to mention that the white spots are due to an upstream issue with Datashader, which will be fixed in the next release.

holoviz/datashader#1329

Can confirm, speckles are gone...though I am getting a lot of:

/Users/mathomp4/installed/Core/GEOSpyD/24.4.0-0_py3.11/2024-06-11/lib/python3.11/site-packages/xarray/namedarray/core.py:496: UserWarning: Duplicate dimension names present: dimensions {'ncontact'} appear more than once in dims=('nf', 'ncontact', 'ncontact'). We do not yet support duplicate dimension names, but we do allow initial construction of the object. We recommend you rename the dims immediately to become distinct, as most xarray functionality is likely to fail silently if you do not. To rename the dimensions you will need to set the ``.dims`` attribute of each variable, ``e.g. var.dims=('x0', 'x1')``.
  warnings.warn(

I can get a picture, so this is working, but is this something we should worry about in future?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Support GEOS Cube Sphere Grids
4 participants