water boundary lines should be linestrings (not unclosed polygons)

[nvkelso]

As reported by Mike Bostock in https://twitter.com/mbostock/status/991715624089874432...

• What did you see? Water layer boundary lines are represented as MultiPolygon.
• What did you expect to see? Water layer boundary lines are represented as MultiLineString, which would render correctly just with strokes (not available for fills) in D3 and related tools.
• What map location are you having problems with? San Francisco coastline.
• Screenshot?
  
  

[nvkelso]

The code that creates the boundary "lines" is here:

vector-datasource/vectordatasource/transform.py

Lines 1300 to 1492 in ed4072a

	def exterior_boundaries(ctx):
	"""
	create new fetures from the boundaries of polygons
	in the base layer, subtracting any sections of the
	boundary which intersect other polygons. this is
	added as a new layer if new_layer_name is not None
	otherwise appended to the base layer.
	the purpose of this is to provide us a shoreline /
	river bank layer from the water layer without having
	any of the shoreline / river bank draw over the top
	of any of the base polygons.
	properties on the lines returned are copied / adapted
	from the existing layer using the new_props dict. see
	_make_new_properties above for the rules.
	buffer_size determines whether any buffering will be
	done to the index polygons. a judiciously small
	amount of buffering can help avoid "dashing" due to
	tolerance in the intersection, but will also create
	small overlaps between lines.
	any features in feature_layers[layer] which aren't
	polygons will be ignored.
	note that the `bounds` kwarg should be filled out
	automatically by tilequeue - it does not have to be
	provided from the config.
	"""
	feature_layers = ctx.feature_layers
	zoom = ctx.nominal_zoom
	base_layer = ctx.params.get('base_layer')
	assert base_layer, 'Missing base_layer parameter'
	new_layer_name = ctx.params.get('new_layer_name')
	prop_transform = ctx.params.get('prop_transform')
	buffer_size = ctx.params.get('buffer_size')
	start_zoom = ctx.params.get('start_zoom', 0)
	snap_tolerance = ctx.params.get('snap_tolerance')
	layer = None
	# don't start processing until the start zoom
	if zoom < start_zoom:
	return layer
	# search through all the layers and extract the one
	# which has the name of the base layer we were given
	# as a parameter.
	layer = _find_layer(feature_layers, base_layer)
	# if we failed to find the base layer then it's
	# possible the user just didn't ask for it, so return
	# an empty result.
	if layer is None:
	return None
	if prop_transform is None:
	prop_transform = {}
	features = layer['features']
	# this exists to enable a dirty hack to try and work
	# around duplicate geometries in the database. this
	# happens when a multipolygon relation can't
	# supersede a member way because the way contains tags
	# which aren't present on the relation. working around
	# this by calling "union" on geometries proved to be
	# too expensive (~3x current), so this hack looks at
	# the way_area of each object, and uses that as a
	# proxy for identity. it's not perfect, but the chance
	# that there are two overlapping polygons of exactly
	# the same size must be pretty small. however, the
	# STRTree we're using as a spatial index doesn't
	# directly support setting attributes on the indexed
	# geometries, so this class exists to carry the area
	# attribute through the index to the point where we
	# want to use it.
	class geom_with_area:
	def __init__(self, geom, area):
	self.geom = geom
	self.area = area
	self._geom = geom._geom
	# STRtree started filtering out empty geoms at some version, so
	# we need to proxy the is_empty property.
	self.is_empty = geom.is_empty
	# create an index so that we can efficiently find the
	# polygons intersecting the 'current' one. Note that
	# we're only interested in intersecting with other
	# polygonal features, and that intersecting with lines
	# can give some unexpected results.
	indexable_features = list()
	indexable_shapes = list()
	for shape, props, fid in features:
	if shape.type in ('Polygon', 'MultiPolygon'):
	# the data comes back clipped from the queries now so we
	# no longer need to clip here
	snapped = shape
	if snap_tolerance is not None:
	snapped = _snap_to_grid(shape, snap_tolerance)
	# geometry collections are returned as None
	if snapped is None:
	continue
	# snapping coordinates and clipping shapes might make the shape
	# invalid, so we need a way to clean them. one simple, but not
	# foolproof, way is to buffer them by 0.
	if not snapped.is_valid:
	snapped = snapped.buffer(0)
	# that still might not have done the trick, so drop any polygons
	# which are still invalid so as not to cause errors later.
	if not snapped.is_valid:
	# TODO: log this as a warning!
	continue
	# skip any geometries that may have become empty
	if snapped.is_empty:
	continue
	indexable_features.append((snapped, props, fid))
	indexable_shapes.append(geom_with_area(snapped, props.get('area')))
	index = STRtree(indexable_shapes)
	new_features = list()
	# loop through all the polygons, taking the boundary
	# of each and subtracting any parts which are within
	# other polygons. what remains (if anything) is the
	# new feature.
	for feature in indexable_features:
	shape, props, fid = feature
	boundary = shape.boundary
	cutting_shapes = index.query(boundary)
	for cutting_item in cutting_shapes:
	cutting_shape = cutting_item.geom
	cutting_area = cutting_item.area
	# dirty hack: this object is probably a
	# superseded way if the ID is positive and
	# the area is the same as the cutting area.
	# using the ID check here prevents the
	# boundary from being duplicated.
	is_superseded_way = \
	cutting_area == props.get('area') and \
	props.get('id') > 0
	if cutting_shape is not shape and \
	not is_superseded_way:
	buf = cutting_shape
	if buffer_size is not None:
	buf = buf.buffer(buffer_size)
	boundary = boundary.difference(buf)
	# filter only linestring-like objects. we don't
	# want any points which might have been created
	# by the intersection.
	boundary = _filter_geom_types(boundary, _LINE_DIMENSION)
	if not boundary.is_empty:
	new_props = _make_new_properties(props, prop_transform)
	new_features.append((boundary, new_props, fid))
	if new_layer_name is None:
	# no new layer requested, instead add new
	# features into the same layer.
	layer['features'].extend(new_features)
	return layer
	else:
	# make a copy of the old layer's information - it
	# shouldn't matter about most of the settings, as
	# post-processing is one of the last operations.
	# but we need to override the name to ensure we get
	# some output.
	new_layer_datum = layer['layer_datum'].copy()
	new_layer_datum['name'] = new_layer_name
	new_layer = layer.copy()
	new_layer['layer_datum'] = new_layer_datum
	new_layer['features'] = new_features
	new_layer['name'] = new_layer_name
	return new_layer

With this bit looking at the dimensions and appending the new features:

vector-datasource/vectordatasource/transform.py

Lines 1463 to 1470 in ed4072a

	# filter only linestring-like objects. we don't
	# want any points which might have been created
	# by the intersection.
	boundary = _filter_geom_types(boundary, _LINE_DIMENSION)
	if not boundary.is_empty:
	new_props = _make_new_properties(props, prop_transform)
	new_features.append((boundary, new_props, fid))

Maybe we need to also recast the result as a line?

[nvkelso]

@zerebubuth Can you take a look, please?