andhuang-CLGX
commented
2 years ago Seems like holoviews / geoviews
Open andhuang-CLGX opened 2 years ago
andhuang-CLGX
commented
2 years ago Seems like holoviews / geoviews
greglucas
commented
2 years ago Are you setting PYPROJ_GLOBAL_CONTEXT=ON? See the note on the new docs about that speed-up: https://scitools.org.uk/cartopy/docs/latest/index.html
However, this looks like a lot of pyproj CRS objects are being created, rather than re-using the same object which would be causing the slowdown. I wonder if we should look at putting a cache on the CRS constructors in addition to the Transformers?
# Many calls like this
for i in range(100):
ax.scatter(..., transform=ccrs.PlateCarree())
# rather than storing the object
pc = ccrs.PlateCarree()
for i in range(100):
ax.plot(..., transform=pc)This is a minimal example:
import geoviews as gv
import cartopy.crs as ccrs
gv.extension("bokeh")
gv.tile_sources.CartoDark() * gv.Points(([-88, -95, -100], [40, 50, 60]), crs=ccrs.PlateCarree())This takes 24 seconds to run on cartopy==0.20, and less than 1 second (0:00:00.183278) on cartopy==0.19.0
leeyupeng119
commented
2 years ago I met the same question.
dopplershift
commented
2 years ago @leeyupeng119 and @andhuang-CLGX: Have you tried setting PYPROJ_GLOBAL_CONTEXT=ON in the environment, like @greglucas suggested, and see if that improves things?
andhuang-CLGX
commented
2 years ago Yep I have tried that. https://github.com/holoviz/geoviews/issues/529#issuecomment-934661990 geoviews may implement their own solution in the near future.
marqh
commented
2 years ago Hi All
this seems like a pretty major performance hit on the 0.19 0.20 transition
does it seem like this is an edge case? is this a performance hit that could be seen (oom) for many cases?
Has any consideration has been given to providing access to the proj transforms for users, enabling them to opt out of the pyproj transforms if they desire?
I'm a bit worried that the documentation note:
The v0.20 release uses pyproj for transformations, which could be slower in some situations. doesn't provide sufficient information given the magnitude of performance hit seen here
I feel that there is more that should be considered here than only PYPROJ_GLOBAL_CONTEXT=ON
many thanks marqh
dopplershift
commented
2 years ago I'm happy to consider those kinds of performance improvements, but any such addition needs to consider the maintainability and sustainability of the project. We have a minimal set of maintainers on this project. Updating to the new PROJ API languished for years until Cartopy's inability to work with current versions of PROJ became such a problem and headache for the community writ large that @snowman2 stepped in to do the update.
So bluntly, without a significant and trusted guarantee that whoever adds direct builds against PROJ will stick around and maintain it long term, I can't see how merging it would be in the best interest of the project. I think the right solution is to figure out where the bottlenecks are and add optimizations (or new APIs) to pyproj.
snowman2
commented
2 years ago Here are performance tips identified in the pyproj development:
PYPROJ_GLOBAL_CONTEXT=ON recommendation exists for cartopy.Transformer class in pyproj.marqh
commented
2 years ago i agree with @dopplershift and the comments on the importance of maintainability
i would also like to recognise the contributions of @snowman2 in delivering the updates to get the code base working with up to date proj
i think that it would be worthwhile exploring potential benefit opportunity, whilst keeping in mind that detailed optimisation tends to lead to complexity which is a maintainability headache.
i wonder what practical steps could be identified to see if this performance can be improved without major code base maintenance issues.
i'd like to help with this endeavour, but i'll state up front that I'm not about to add direct builds against PROJ and commit to long term maintenance of said builds myself
many thanks all marqh
greglucas
commented
2 years ago @marqh, investigations into the slowdown would be most welcome!
One of the biggest areas that could use a speed-up would be the bisection routine. https://github.com/SciTools/cartopy/blob/72216fd3d68aaf3fb693327ff88baa032a0170db/lib/cartopy/trace.pyx#L441 Calling pyproj's Transform with a single point over and over is not ideal. So, thinking of a way to avoid that iteration, or speed it up and call the transform with 4-10 points each time as a look-ahead, may be worthwhile.
For caching the CRS constructor, I am not sure whether that should be implemented on Cartopy's side, or if a PR to pyproj would be the better location for it. You'd likely need to override __new__ to return from a global cache stored somewhere. It seems doable if the same parameters are passed in to the constructor, but that could be optimistic thinking too.
There are already performance metrics that can be run in the benchmarks folder: https://github.com/SciTools/cartopy/tree/master/benchmarks
those can be run with asv. We aren't currently running them on a CI system, but perhaps something could be set up to check for major regressions with it automated.
leeyupeng119
commented
2 years ago @leeyupeng119 and @andhuang-CLGX: Have you tried setting
PYPROJ_GLOBAL_CONTEXT=ONin the environment, like @greglucas suggested, and see if that improves things?
Yes, I have try, but still slow.
akrherz
commented
2 years ago I'm digging at this general issue too now and find that PYPROJ_GLOBAL_CONTEXT=ON is not making any difference in the run-time. For an interesting note, my current debugging script goes from a run-time of 4 seconds to 11 seconds just by adding cfeature.LAND to the plot. That's a bit puzzling and where I am currently looking.
ktyle
commented
2 years ago @akrherz I find an even larger performance penalty by adding cfeature.OCEAN (esp. with 10m res and transforming to another CRS), but this pre-dates 0.20.
akrherz
commented
2 years ago @ktyle Yeah, I realized that the ocean is not a necessary feature to draw. Just make the axes background the ocean color and then draw the land, which is much less expensive to do.
So I think I just reproduced what @greglucas stated above, the segment-by-segment reprojection in trace.pyx is what is very slow and what causes cfeature.LAND to be so slow, given the large number of segments it has.
greglucas
commented
2 years ago @akrherz, can you see if PR https://github.com/SciTools/cartopy/pull/1918 helps your use-case at all?
akrherz
commented
2 years ago @akrherz, can you see if PR #1918 helps your use-case at all?
Sorry for the delay, sadly not much.
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(111, projection=ccrs.Mercator())
feature = cfeature.NaturalEarthFeature(
name="ocean",
category="physical",
scale="10m",
edgecolor="#000000",
facecolor="#AAAAAA",
)
ax.add_feature(feature)
fig.savefig('test.png')~13.8 seconds to ~13.4 seconds.
maximlt
commented
2 years ago Just to report the output of a tiny benchmark, comparing a projection creation and a projection comparison with a simple string that is expected to return False:
0.19.0.post1:
0.20.1:
In this case, the creation of a projection is about 25x slower and the equality test is 3,500,000x slower.
dopplershift
commented
2 years ago It's unclear from your example if there's something else going on. You run p = cartopy.crs.PlateCarree(), but your timeit is checking proj.
maximlt
commented
2 years ago @dopplershift indeed this might have been misleading. I updated the screenshots in my previous message. The outcome is still the same.
snowman2
commented
2 years ago @maximlt, the updated method of PROJ is more robust when comparing CRS and has many different methods to create a CRS (PROJ string, Authority, WKT, PROJ JSON, etc...). Since proj is not a valid input to create a CRS, PROJ will exhaust all possibilities to create a CRS for you every time you perform a comparison in the timeit operation. I would imagine if you were able to successfully create a cartopy.crs.CRS beforehand for comparison purposes, it would be more efficient.
snowman2
commented
2 years ago 
dopplershift
commented
2 years ago While the microbenchmark definitely agrees with the hotspots identified by @andhuang-CLGX's profiling, I think more detail at what holoviews is doing is necessary to see how equals and __init__ to know exactly what's causing it to be slow.
maximlt
commented
2 years ago Thanks for your replies and explanations! I'm well aware that using pyproj is an improvement and reported this small benchmark here for others to be aware of that, since geoviews (actually holoviews) might not be the only library impacted by this change.
warrickball
commented
2 years ago I was going to post a new issue but there's reasonable chance mine is related to this.
I'm trying to place markers on a simple world map and am finding the map takes a very long time (about 30s) to zoom in on, say, Great Britain, using Matplotlib's standard zoom tool (the magnifying glass icon). This happens in both the Robinson or Mercator projections. Here's the plot where I'm trying to zoom in:
import matplotlib.pyplot as pl
import cartopy
ax = pl.subplot(1, 1, 1, projection=cartopy.crs.Robinson())
ax.set_global()
ax.add_feature(cartopy.feature.OCEAN, zorder=0)
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
pl.show()I'm running Fedora 35 and installed geos 3.9.2 and proj 8.2.0 through the default repos. I have Python 3.10.1 and installed cartopy with python3 -m pip install cartopy --user, which has given me version 0.20.1. From my matplotlibrc, I think I'm using the TkAgg backend for matplotlib, which is version 3.5.1.
I'll keep experimenting to see if anything jumps out at me. This might be related to #1895 but PYPROJ_GLOBAL_CONTEXT=ON didn't help and I got build errors when downgrading to cartopy<0.20.
warrickball
commented
2 years ago I tried running my script through cProfile by opening the plot, zooming and then pressing the close button, so that the program exited as soon as the zoom (on a region of northern Europe, including e.g. North Sea and Baltic Sea) was finished. Not sure if it might be helpful. The full terminal output is attached here and the snippet below is the first few lines (until tottime < 1s).
209131486 function calls (208964453 primitive calls) in 170.533 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
11193149 27.966 0.000 30.714 0.000 coords.py:61(__iter__)
9112 25.372 0.003 28.662 0.003 linestring.py:229(geos_linestring_from_py)
49192 9.034 0.000 11.600 0.000 coords.py:164(xy)
3583410 8.797 0.000 15.593 0.000 {method '_transform' of 'pyproj._transformer._Transformer' objects}
7166820 8.612 0.000 17.505 0.000 utils.py:88(_copytobuffer)
63741 7.517 0.000 25.296 0.000 polygon.py:438(geos_linearring_from_py)
3583410 6.864 0.000 44.053 0.000 transformer.py:649(transform)
1017607 6.083 0.000 19.258 0.000 coords.py:76(__getitem__)
14552 6.048 0.000 50.662 0.003 {cartopy.trace.project_linear}
7166636 4.858 0.000 4.858 0.000 utils.py:66(_copytobuffer_return_scalar)
33325097 4.201 0.000 4.201 0.000 {built-in method _ctypes.byref}
37087951 4.061 0.000 4.063 0.000 {built-in method builtins.isinstance}
2184301 3.451 0.000 11.461 0.000 coords.py:43(_update)
2350772 3.287 0.000 5.213 0.000 predicates.py:23(__call__)
7166820 2.539 0.000 2.539 0.000 utils.py:138(_convertback)
1 2.522 2.522 168.790 168.790 {built-in method exec}
2350715 2.464 0.000 8.470 0.000 base.py:709(is_empty)
19 2.256 0.119 62.185 3.273 crs.py:955(_attach_lines_to_boundary)
11078928 2.145 0.000 3.286 0.000 linestring.py:258(_coords)
1116973 2.047 0.000 7.965 0.000 coords.py:51(__len__)
3609193 1.916 0.000 2.628 0.000 enum.py:359(__call__)
3583413 1.907 0.000 4.507 0.000 enums.py:13(create)
14394652 1.675 0.000 3.156 0.000 {built-in method builtins.hasattr}
3583412 1.552 0.000 1.552 0.000 transformer.py:313(_transformer)
250 1.490 0.006 1.490 0.006 {method 'read' of '_ssl._SSLSocket' objects}
9464100 1.274 0.000 1.274 0.000 base.py:228(_geom)
1486562/1465361 1.105 0.000 1.123 0.000 base.py:245(__setattr__)
6494596 1.028 0.000 1.028 0.000 {method 'append' of 'array.array' objects}
...This shapely issue might be related to my issue, though it also makes me suspect that my issue isn't related to the original one in this thread. The runtime of my example decreased to ~12s when I installed shapely from the GitHub repo following the example in this comment:
pip install --verbose --no-binary=shapely "git+https://github.com/shapely/shapely@maint-1.8#egg=shapely"Happy to open a new issue if mine is independent.
jaredalee
commented
2 years ago I realized that the ocean is not a necessary feature to draw. Just make the axes background the ocean color and then draw the land, which is much less expensive to do. @akrherz
I've also had numerous problems with Cartopy 0.20.1/Shapely 1.8 suddenly being unbearably slow after I recently upgraded to them. Plotting the terrain field for a large WRF domain of 1062 x 512 points (dx = 1 km) suddenly took nearly 11m 30s. That's unacceptably slow. With my prior version of Cartopy (0.16, I think??), it took "only" 7 min, which was still far too slow to be useful. (I'm running this on NCAR's Cheyenne, on a head node.)
I set PYPROJ_GLOBAL_CONTEXT=ON, but that had little effect on the speed for this plot. (I noticed more of a speed-up for other plots I tried of different data, though.)
I then used the transform_first argument to plt.contourf ([https://scitools.org.uk/cartopy/docs/latest/gallery/scalar_data/contour_transforms.html]), and that resulted in a noticeable speed-up, to 10m 45s.
However, commenting out the call to draw the oceans (ax.add_feature(oceans)), but instead assigning my ocean color using ax.set_facecolor, resulted in a SIGNIFICANT speed-up, all the way down to 12 s. If I eliminated the transform_first argument again, it slowed down to just over 1 min, so transform_first is still quite helpful.
So for a temporary workaround to anyone experiencing similar major slowdowns with Cartopy 0.20+, hopefully implementing these steps will help:
transform_first keyword in the contourf call (moderate speed-up).PYPROJ_GLOBAL_CONTEXT=ON for a potential speed-up (in some situations, but not all).I also very much look forward to Cartopy 0.20.2 release to eliminate the ShapelyDeprecationErrors (#1936) that also started happening constantly once I upgraded.
warrickball
commented
2 years ago However, commenting out the call to draw the oceans (ax.add_feature(oceans)), but instead assigning my ocean color using ax.set_facecolor, resulted in a SIGNIFICANT speed-up, all the way down to 12 s.
Thanks @jaredalee! I found this very helpful, reducing the time to zoom in on my point map from ~30s to about ~3s.
If anyone else is curious, Cartopy's water colour is np.array((152, 183, 226)) / 256.:
greglucas
commented
2 years ago Thanks for posting these results! I think this points to the need for some documentation around potential speedups that people can try all in one place rather than spread throughout the GitHub Issues.
From my reading of everything, it seems like the major issue is the oceans feature and not land or coastlines? This may mean we should look into optimizing polygon edge closing with the boundaries.
You should be able to get the colors programmatically from that location too: cartopy.feature.COLORS["water"] https://scitools.org.uk/cartopy/docs/latest/reference/generated/cartopy.feature.COLORS.html#cartopy.feature.COLORS
warrickball
commented
2 years ago D'oh! Thanks for the better programmatic access.
While this is definitely faster, it still gets slow as I zoom in further, even without the ocean. Here's my current example:
#!/usr/bin/env python
import matplotlib.pyplot as pl
import cartopy
ax = pl.subplot(1, 1, 1, projection=cartopy.crs.Robinson())
ax.set_global()
ax.set_facecolor(cartopy.feature.COLORS['water'])
ax.add_feature(cartopy.feature.LAND, zorder=0, edgecolor='black')
pl.show()
It took about a minute to zoom in on a region roughly corresponding to Scotland (including the Hebrides and a bit of Northern Ireland). Here's the start of the cProfile output:
78692232 function calls (78587675 primitive calls) in 62.467 seconds
Ordered by: internal time
ncalls tottime percall cumtime percall filename:lineno(function)
1 5.294 5.294 60.703 60.703 {built-in method exec}
27312 4.976 0.000 6.410 0.000 coords.py:164(xy)
1703573 4.340 0.000 7.629 0.000 {method '_transform' of 'pyproj._transformer._Transformer' objects}
3407146 4.220 0.000 8.479 0.000 utils.py:88(_copytobuffer)
35554 3.796 0.000 12.785 0.000 polygon.py:438(geos_linearring_from_py)
1703573 3.309 0.000 21.418 0.000 transformer.py:649(transform)
493202 3.052 0.000 9.624 0.000 coords.py:76(__getitem__)
6830 3.008 0.000 24.696 0.004 {cartopy.trace.project_linear}
3406556 2.308 0.000 2.308 0.000 utils.py:66(_copytobuffer_return_scalar)
1068583 1.767 0.000 5.798 0.000 coords.py:43(_update)
1197644 1.728 0.000 2.754 0.000 predicates.py:23(__call__)
12855265 1.471 0.000 1.475 0.000 {built-in method builtins.isinstance}
1197596 1.288 0.000 4.454 0.000 base.py:709(is_empty)
3407146 1.259 0.000 1.259 0.000 utils.py:138(_convertback)
413541 1.108 0.000 1.212 0.000 coords.py:61(__iter__)
548018 1.041 0.000 4.053 0.000 coords.py:51(__len__)
1729426 0.918 0.000 1.252 0.000 enum.py:359(__call__)
1703576 0.883 0.000 2.107 0.000 enums.py:13(create)
6429967 0.866 0.000 0.866 0.000 {built-in method _ctypes.byref}
6856486 0.813 0.000 1.395 0.000 {built-in method builtins.hasattr}
1703575 0.742 0.000 0.742 0.000 transformer.py:313(_transformer)
...I'm not sure if I should expect Cartopy to take that long but it surprises me.
ktyle
commented
2 years ago Thanks @jaredalee and @warrickball ! @greglucas , including the Oceans feature has always involved a big performance hit, even pre-dating 0.20 in my experience. Using the ax.set_facecolor workaround makes a huge difference.
I am still investigating but I found that somewhere in my app, the projecting is taking 97% of the load time.
Will update this as I figure out where it's coming from.
Some context, I am using geoviews, datashader to plot millions of points. I am not sure if geoviews / datashader is using an outdated cartopy method to transform the points; will test.