geoviews.feature extents

[scaine1]

It would be nice if the extents given to the geoviews.feature determine the extent of data that is extracted/used, rather than just determing the viewable window size.

Below is an example that plots some fake data over a small region. Here the window extent is chosen by a) the dataset used and b) the extents given to the geoviews.feature.coastline

However, if you zoom out you can see that the coastline data is actually available for the entire globe rather than just the subdomain I have chosen.

This has implications for the size of files generated (HoloMaps) or data downloaded (DynamicMaps). This is especially noticiable when using high resolution data. If you turn on the use of 10m coastline resolution, you will see that the file size goes from ~approximately 150k to 8.8M. It is my suspiciion that the majority of the increased data comes from the highresolution data in regions on the world we don't actually care about.

import numpy as np
import xarray as xr
import holoviews as hv
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs as ccrs

def create_data():
    """
    create fake example data
    """
    times = [np.datetime64('2018-02-22 00:00:00')]
    lons = np.arange(144, 150, 1)
    lats = np.arange(-45, -39, 1)

    dims = ('time', 'latitude', 'longitude')
    shape = (len(times), len(lats), len(lons))
    coords = {'latitude':lats, 'longitude':lons, 'time':times}

    pres = np.zeros(shape)
    pres[0, :, :] = 10
    pres_xr = xr.DataArray(pres, dims=dims, coords=coords)
    xr_ds = xr.Dataset({'time': times, 'latitude': lats,
        'longitude': lons,  'pres': (dims, pres)})
    return xr_ds

hv.extension('bokeh')
renderer = hv.renderer('bokeh')
hv.output('size=200')
hv.opts("Image (cmap='Spectral_r')")

xr_ds = create_data()
vdims = ['pres']
kdims= ['time', 'latitude', 'longitude']
gv_ds = gv.Dataset(xr_ds, kdims=kdims, vdims=vdims, crs=ccrs.PlateCarree())

hv.Dimension.type_formatters[np.datetime64] = '%Y-%m-%d %H'

geo_dims = ['longitude', 'latitude']

extents = (144.0, -45.0, 149.0, -40.0)
# using extents on the geo feature sets the plot window extents,
# however, all the data is still extracted/used

#hv.opts("Feature [scale='10m']")
coast = gf.coastline.clone(extents=extents)
img = gv_ds.to(gv.Image, geo_dims) * coast
renderer.save(img, 'test_gf')

[tsch0nn]

I will have use for this enhancement as well!

[ahuang11]

This seems to work decently, but I don't know where to dump it. I was looking at FeaturePlot, get_extents(), but I don't know how to get the extents set.

import geoviews as gv
from shapely.geometry import box, MultiLineString
from shapely.ops import split, unary_union

bbox = box(-180, 0, 180, 90)

x = gv.feature.coastline.geoms().geom()
geoms = [geom.intersection(bbox) for geom in x]
unary_union(geoms)

[poplarShift]

I think GeoViews should (but doesn't at the moment) have a way of extracting Feature data at different resolutions, setting it directly at the data level instead of doing it through the scale option, which is what seems to be the problem here, if I understand correctly. But just getting the data directly from the shapefile gives me a small html file of ~63 kb, @ahuang11:

import numpy as np
import xarray as xr
import holoviews as hv
from holoviews import opts
import geoviews as gv
import geoviews.feature as gf
from cartopy import crs as ccrs
import cartopy.io.shapereader as shpreader

from shapely.geometry import box, MultiLineString
from shapely.ops import split, unary_union

def create_data():
    """
    create fake example data
    """
    times = [np.datetime64('2018-02-22 00:00:00')]
    lons = np.arange(144, 150, 1)
    lats = np.arange(-45, -39, 1)

    dims = ('time', 'latitude', 'longitude')
    shape = (len(times), len(lats), len(lons))
    coords = {'latitude':lats, 'longitude':lons, 'time':times}

    pres = np.zeros(shape)
    pres[0, :, :] = 10
    pres_xr = xr.DataArray(pres, dims=dims, coords=coords)
    xr_ds = xr.Dataset({'time': times, 'latitude': lats,
        'longitude': lons,  'pres': (dims, pres)})
    return xr_ds

xr_ds = create_data()
vdims = ['pres']
kdims= ['time', 'latitude', 'longitude']
gv_ds = gv.Dataset(xr_ds, kdims=kdims, vdims=vdims, crs=ccrs.PlateCarree())

hv.Dimension.type_formatters[np.datetime64] = '%Y-%m-%d %H'

geo_dims = ['longitude', 'latitude']

extents = (144.0, -45.0, 149.0, -40.0)

renderer = hv.renderer('bokeh')
hv.output('size=200')
hv.opts("Image (cmap='Spectral_r')")

bbox = box(*extents)

shpfilename = shpreader.natural_earth(resolution='10m',
                                      category='physical',
                                      name='coastline')

reader = shpreader.Reader(shpfilename)
coastline = unary_union([geom.geometry.intersection(bbox) for geom in reader.records()])

img = gv_ds.to(gv.Image, geo_dims) * gv.Shape(coastline)
renderer.save(img, 'test_gf')

[philippjfr]

I think GeoViews should (but doesn't at the moment) have a way of extracting Feature data at different resolutions, setting it directly at the data level instead of doing it through the scale option, which is what seems to be the problem here

That's not really accurate, you can always access the Feature element's data like this and request a specific scale using the with_scale method, putting it all together this is the easiest approach:

from shapely.ops import unary_union

extents = (144.0, 149.0, -45.0, -40.0)
coastline = gv.feature.coastline.data.with_scale('10m')
unary_union(list(coastline.intersecting_geometries(extents)))

[poplarShift]

Nice! I wasn't aware of .with_scale. But you'd still pass that into Shape instead of setting the extent directly on the Feature, if I'm not mistaken.

[philippjfr]

Yes, we definitely still want to expose this functionality directly on Feature in some way, haven't yet decided what the best API would be. A plot option would be easy but this could also work via the extents attribute on the element or even via .select and __getitem__.

[poplarShift]

What about setting it through a Dimension(range=...) kind of declaration on the Feature element? Too convoluted?

[philippjfr]

Setting the dimension range does not have the right semantics, it defines the extent of the domain but generally does not slice the actual data. I'm generally in favor of eventually removing the extents parameter on elements but this is one case where I think it at least has useful semantics which aren't already covered by dimension ranges and the xlim/ylim options.

[zxdawn]

@poplarShift @philippjfr I have the similar problem when adding shpfile. This is part of the script:

shpfile = '/public/home/zhangxin/models/data/wrfpython/data/CHN_adm_shp/CHN_adm3.shp'
extents = (122.5, 39.5, 126.5, 44.5)
bbox = box(*extents)
reader = shpreader.Reader(shpfile)
coastline = unary_union([geom.geometry.intersection(bbox) for geom in reader.records()])
.......
gv_data = gv_data.opts(cmap=cmap, colorbar=True, width=600, height=500, xlabel='Longitude', ylabel='Latitude')
gv_data = (gv_data * gv.Shape(coastline)).opts(opts.Feature(fill_color='none', color='white', line_color='white'))
gv.save(gv_data, save_path + 'chn_FY4A_Cha.png', fmt='png')

This is the result:

There're three problems:

1. How to show the shp lines without filled color? The left and right vertical line shouldn't be shown. I need to remove it as shp file doesn't include it.
2. How to remove icons on the right corner?
3. The yaxis label name is wrong ...

[ahuang11]

import geopandas as gpd
import geoviews as gv
gv.extension('bokeh')

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

gv.Polygons(world).opts('Polygons', toolbar='disable', fill_alpha=0)

1. fill_alpha=0 removes the fill, but not sure what you mean by the left and right vertical lines because that's just the edge color / line?

2. toolbar='disable' or use matplotlib backend

3. Not sure about 3; could that be a locale thing, i.e. if you're using a Chinese computer, so it supports writing Chinese characters vertically like that, which somehow affects English letters as well?

[zxdawn]

@ahuang11 Thanks! How about just plotting part of coastlines based on the extent of data? Plotting the whole shp costs a lot of time.

[ahuang11]

See above I think.

[zxdawn]

@ahuang11 OK. This one works:

shpfile = '/public/home/zhangxin/models/data/wrfpython/data/CHN_adm_shp/CHN_adm3.shp'
extents = (122.5, 39.5, 126.5, 44.5)
bbox = box(*extents)
reader = shpreader.Reader(shpfile)
coastline = gv.Shape(unary_union([geom.geometry.intersection(bbox) for geom in reader.records()]))
shp = gv.Polygons(coastline).opts('Polygons', toolbar='disable', fill_alpha=0)
gv_data = gv_data * shp

You can see the issue of these two vertical lines caused by extents. That is more clear when the extents is changed to extents = (125, 35, 130, 45):

[ahuang11]

Sorry, I'm still having trouble understanding. Can you circle the parts you're looking at?

This line?

[zxdawn]

@ahuang11 Yes, that line. I guess it close the shape automatically by extent and coastline.

[ahuang11]

You could set a tighter xlim/ylim (zoom) while a coarser bbox (crop) so you don't see that vertical line and the size is lessened.

[ahuang11]

http://geopandas.org/indexing.html I just discovered there's a.cx method which might be useful here.

[ahuang11]

I would like to help make a PR for this, but I don't know where/how to start.

I was thinking of adding a subset or crop kwarg like this: gv.feature.coastline(subset=(x1, x2, y1, y2)) sort of similar to the scale kwarg gv.feature.coastline(scale='50m') so I tried searching for the scale kwarg in the code, but I only found it in https://github.com/holoviz/geoviews/blob/master/geoviews/element/geo.py#L152 which is only called when doing gv.feature.coastline(scale='50m').geoms

Any ideas appreciated!