Improve x-limit, y-limit, and rasterize support for geo features (#1135)

examples/reference/xarray/vectorfield.ipynb

@@ -24,7 +24,6 @@
    "source": [
     "import numpy as np\n",
     "import xarray as xr\n",
-    "import geoviews as gv\n",
     "import cartopy.crs as ccrs"
    ]
   },
@@ -88,7 +87,24 @@
    "outputs": [],
    "source": [
     "ds.hvplot.vectorfield(x='x', y='y', angle='angle', mag='mag',\n",
-    "                      hover=False, geo=True).opts(magnitude='mag') * gv.tile_sources.CartoLight"
+    "                      hover=False, geo=True, tiles=\"CartoLight\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If you set `coastline` or `features` it will keep the original crs and transform the features to the data crs."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ds.hvplot.vectorfield(x='x', y='y', angle='angle', mag='mag',\n",
+    "                      hover=False, geo=True, coastline=True)"
    ]
   }
  ],

hvplot/converter.py

@@ -446,20 +446,19 @@ def __init__(
                         "Projection must be defined as cartopy CRS or "
                         f"one of the following CRS string:\n {all_crs}")
 
-            projection = projection or (ccrs.GOOGLE_MERCATOR if tiles else self.crs)
-            if tiles and projection != ccrs.GOOGLE_MERCATOR:
+            self.output_projection = projection or (ccrs.GOOGLE_MERCATOR if tiles else self.crs)
+            if tiles and self.output_projection != ccrs.GOOGLE_MERCATOR:
                 raise ValueError(
                     "Tiles can only be used with output projection of "
                     "`cartopy.crs.GOOGLE_MERCATOR`. To get rid of this error "
                     "remove `projection=` or `tiles=`"
                 )
-
-            if self.crs != projection:
+            if self.crs != projection and (xlim or ylim):
                 px0, py0, px1, py1 = ccrs.GOOGLE_MERCATOR.boundary.bounds
                 x0, x1 = xlim or (px0, px1)
                 y0, y1 = ylim or (py0, py1)
                 extents = (x0, y0, x1, y1)
-                x0, y0, x1, y1 = project_extents(extents, self.crs, projection)
+                x0, y0, x1, y1 = project_extents(extents, self.crs, self.output_projection)
                 if xlim:
                     xlim = (x0, x1)
                 if ylim:
@@ -1363,6 +1362,10 @@ def method_wrapper(ds, x, y):
             if self._dim_ranges.get('c', (None, None)) != (None, None):
                 style['clim'] = self._dim_ranges['c']
 
+        if self.geo and self.crs != self.output_projection:
+            import geoviews as gv
+            obj = gv.project(obj, projection=self.output_projection)
+
         processed = operation(obj, **opts)
 
         if self.dynspread:
@@ -1389,7 +1392,7 @@ def _apply_layers(self, obj):
                 param.main.param.warning(
                     "coastline scale of %s not recognized, must be one "
                     "'10m', '50m' or '110m'." % self.coastline)
-            obj = obj * coastline
+            obj = obj * coastline.opts(projection=self.output_projection)
 
         if self.features:
             import geoviews as gv
@@ -1411,32 +1414,46 @@ def _apply_layers(self, obj):
                     else:
                         feature_obj = feature_obj.opts(scale=scale)
                 if feature_obj.group in ["Land", "Ocean"]:
-                    obj = feature_obj * obj  # Underlay land/ocean
+                    # Underlay land/ocean
+                    obj = feature_obj.opts(projection=self.output_projection) * obj
                 else:
-                    obj = obj * feature_obj  # overlay everything else
-
-        if self.tiles:
-            tile_source = 'EsriImagery' if self.tiles == 'ESRI' else self.tiles
-            warning = ("{} tiles not recognized, must be one of: {} or a tile object".format(tile_source, sorted(hv.element.tile_sources)))
-            if tile_source is True:
-                tiles = hv.element.tiles.OSM()
-            elif tile_source in hv.element.tile_sources.keys():
-                tiles = hv.element.tile_sources[tile_source]()
-            elif tile_source in hv.element.tile_sources.values():
-                tiles = tile_source()
-            elif isinstance(tile_source, hv.element.tiles.Tiles):
-                tiles = tile_source
-            elif self.geo:
-                from geoviews.element import WMTS
-                if isinstance(tile_source, WMTS):
-                    tiles = tile_source
-                else:
-                    param.main.param.warning(warning)
-            else:
-                param.main.param.warning(warning)
+                    # overlay everything else
+                    obj = obj * feature_obj.opts(projection=self.output_projection)
+
+        if self.tiles and not self.geo:
+            tiles = self._get_tiles(
+                self.tiles,
+                hv.element.tile_sources,
+                hv.element.tiles.Tiles
+            )
+            obj = tiles * obj
+        elif self.tiles and self.geo:
+            import geoviews as gv
+            tiles = self._get_tiles(
+                self.tiles,
+                gv.tile_sources.tile_sources,
+                (gv.element.WMTS, hv.element.tiles.Tiles),
+            )
             obj = tiles * obj
         return obj
 
+    def _get_tiles(self, source, sources, types):
+        tile_source = 'EsriImagery' if self.tiles == 'ESRI' else self.tiles
+        if tile_source is True:
+            tiles = sources["OSM"]()
+        elif tile_source in sources:
+            tiles = sources[tile_source]()
+        elif tile_source in sources.values():
+            tiles = tile_source()
+        elif isinstance(tile_source, types):
+            tiles = tile_source
+        else:
+            msg = (
+                f"{tile_source} tiles not recognized, must be one of: {sorted(sources)} or a tile object"
+            )
+            raise ValueError(msg)
+        return tiles
+
     def _merge_redim(self, ranges, attr='range'):
         redim = dict(self._redim)
         for k, r in ranges.items():

hvplot/tests/testgeo.py

@@ -111,6 +111,43 @@ def test_plot_with_projection_raises_an_error_when_tiles_set(self):
         with self.assertRaisesRegex(ValueError, "Tiles can only be used with output projection"):
             da.hvplot.image('x', 'y', crs=self.crs, projection='Robinson', tiles=True)
 
+    def test_overlay_with_projection(self):
+        # Regression test for https://github.com/holoviz/hvplot/issues/1090
+        df = pd.DataFrame({"lon": [0, 10], "lat": [40, 50], "v": [0, 1]})
+
+        plot1 = df.hvplot.points(x="lon", y="lat", s=200, c="y", geo=True, tiles="CartoLight")
+        plot2 = df.hvplot.points(x="lon", y="lat", c="v", geo=True)
+
+        # This should work without erroring
+        plot = plot1 * plot2
+        hv.renderer("bokeh").get_plot(plot)
+
+    def test_geo_with_rasterize(self):
+        import xarray as xr
+        import cartopy.crs as ccrs
+        import geoviews as gv
+        try:
+            from holoviews.operation.datashader import rasterize
+        except:
+            raise SkipTest('datashader not available')
+
+        ds = xr.tutorial.open_dataset("air_temperature")
+        hvplot_output = ds.isel(time=0).hvplot.points(
+            "lon",
+            "lat",
+            crs=ccrs.PlateCarree(),
+            projection=ccrs.LambertConformal(),
+            rasterize=True,
+            dynamic=False,
+            aggregator="max",
+        )
+
+        p1 = gv.Points(ds.isel(time=0), kdims=["lon", "lat"], crs=ccrs.PlateCarree())
+        p2 = gv.project(p1, projection=ccrs.LambertConformal())
+        expected = rasterize(p2, dynamic=False, aggregator="max")
+
+        xr.testing.assert_allclose(hvplot_output.data, expected.data)
+
 
 class TestGeoAnnotation(TestCase):
 
@@ -141,27 +178,55 @@ def test_plot_with_coastline_sets_geo_by_default(self):
     def test_plot_with_coastline_scale(self):
         plot = self.df.hvplot.points('x', 'y', geo=True, coastline='10m')
         opts = plot.get(1).opts.get('plot')
-        self.assertEqual(opts.kwargs, {'scale': '10m'})
+        assert opts.kwargs["scale"] == '10m'
 
     def test_plot_with_tiles(self):
-        plot = self.df.hvplot.points('x', 'y', geo=True, tiles=True)
+        plot = self.df.hvplot.points('x', 'y', geo=False, tiles=True)
         self.assertEqual(len(plot), 2)
         self.assertIsInstance(plot.get(0), hv.Tiles)
         self.assertIn('openstreetmap', plot.get(0).data)
 
+    def test_plot_with_tiles_with_geo(self):
+        import geoviews as gv
+
+        plot = self.df.hvplot.points('x', 'y', geo=True, tiles=True)
+        self.assertEqual(len(plot), 2)
+        self.assertIsInstance(plot.get(0), gv.element.WMTS)
+        self.assertIn('openstreetmap', plot.get(0).data)
+
     def test_plot_with_specific_tiles(self):
-        plot = self.df.hvplot.points('x', 'y', geo=True, tiles='ESRI')
+        plot = self.df.hvplot.points('x', 'y', geo=False, tiles='ESRI')
         self.assertEqual(len(plot), 2)
         self.assertIsInstance(plot.get(0), hv.Tiles)
         self.assertIn('ArcGIS', plot.get(0).data)
 
+    def test_plot_with_specific_tiles_geo(self):
+        import geoviews as gv
+        plot = self.df.hvplot.points('x', 'y', geo=True, tiles='ESRI')
+        self.assertEqual(len(plot), 2)
+        self.assertIsInstance(plot.get(0), gv.element.WMTS)
+        self.assertIn('ArcGIS', plot.get(0).data)
+
     def test_plot_with_specific_tile_class(self):
-        plot = self.df.hvplot.points('x', 'y', geo=True, tiles=hv.element.tiles.EsriImagery)
+        plot = self.df.hvplot.points('x', 'y', geo=False, tiles=hv.element.tiles.EsriImagery)
         self.assertEqual(len(plot), 2)
         self.assertIsInstance(plot.get(0), hv.Tiles)
         self.assertIn('ArcGIS', plot.get(0).data)
 
+    def test_plot_with_specific_tile_class_with_geo(self):
+        import geoviews as gv
+        plot = self.df.hvplot.points('x', 'y', geo=True, tiles=gv.tile_sources.EsriImagery)
+        self.assertEqual(len(plot), 2)
+        self.assertIsInstance(plot.get(0), gv.element.WMTS)
+        self.assertIn('ArcGIS', plot.get(0).data)
+
     def test_plot_with_specific_tile_obj(self):
+        plot = self.df.hvplot.points('x', 'y', geo=False, tiles=hv.element.tiles.EsriImagery())
+        self.assertEqual(len(plot), 2)
+        self.assertIsInstance(plot.get(0), hv.Tiles)
+        self.assertIn('ArcGIS', plot.get(0).data)
+
+    def test_plot_with_specific_tile_obj_with_geo(self):
         plot = self.df.hvplot.points('x', 'y', geo=True, tiles=hv.element.tiles.EsriImagery())
         self.assertEqual(len(plot), 2)
         self.assertIsInstance(plot.get(0), hv.Tiles)
@@ -294,10 +359,16 @@ def test_points_hover_cols_with_by_set_to_name(self):
             assert element.vdims == []
 
     def test_points_project_xlim_and_ylim(self):
+        points = self.cities.hvplot(geo=False, xlim=(-10, 10), ylim=(-20, -10))
+        opts = hv.Store.lookup_options('bokeh', points, 'plot').options
+        np.testing.assert_equal(opts['xlim'], (-10, 10))
+        np.testing.assert_equal(opts['ylim'], (-20, -10))
+
+    def test_points_project_xlim_and_ylim_with_geo(self):
         points = self.cities.hvplot(geo=True, xlim=(-10, 10), ylim=(-20, -10))
         opts = hv.Store.lookup_options('bokeh', points, 'plot').options
-        assert opts['xlim'] == (-10, 10)
-        assert opts['ylim'] == (-20, -10)
+        np.testing.assert_allclose(opts['xlim'], (-10, 10))
+        np.testing.assert_allclose(opts['ylim'], (-20, -10))
 
     def test_polygons_by_subplots(self):
         polygons = self.polygons.hvplot(geo=True, by="name", subplots=True)