nst-guide / web

Mapping website for National Scenic Trails
https://nst.guide
GNU General Public License v3.0
5 stars 0 forks source link

PCT Elevation Chart #16

Open kylebarron opened 4 years ago

kylebarron commented 4 years ago

Updating elevation chart:

Tippecanoe only stores two dimensions when generating vector tiles. So if I take the 3d lines and convert them to vector tiles, I'll have no elevation data.

So my options are essentially:

  1. GeoJSON with 3 dimensions
  2. vector tiles of points with metadata on each point
  3. GeoJSON of points with metadata on each

I think 2 is the best, because I can also encode the trail mile in the data.

Note that a drawback of this is that it could simplify the geometry when zoomed out, and thus you wouldn't have an exact elevation chart at low zooms. You could only show at a higher zoom, and also there's probably a tippecanoe option to not simplify as much.

Given this drawback, alternatively, 3 might be the best because there would be no simplification. Then you could have the trail line be vector tile

Also to think about:

kylebarron commented 4 years ago

Note that if you run map.queryRenderedFeatures with no geometry provided, by default it runs it for the current viewport. You can still pass a specific layer id, so that you only find elevation features.

kylebarron commented 4 years ago

Code to generate geojson of points with elevation:

import geojson
import numpy as np
from shapely.geometry import LineString, Point
from shapely.ops import linemerge

import geom
from data_source import Halfmile

hm = Halfmile()
gdf = hm.trail_full(False)
l = linemerge(gdf.geometry.values)
l_repr = geom.reproject(l, geom.WGS84, geom.CA_ALBERS)
multiplier = 4282464 / l_repr.length

marginal_dists = []
for ind in range(len(l_repr.coords)):
    if ind == 0:
        marginal_dists.append(0.0)
        continue

    coord_next = l_repr.coords[ind]
    coord_prev = l_repr.coords[ind - 1]
    line = LineString([Point(coord_next), Point(coord_prev)])
    dist = line.length
    marginal_dists.append(dist)

dists_new = np.cumsum(marginal_dists)
# Adjust to actual PCT length
dists_new *= multiplier

features = []
for coord, dist in zip(l.coords, dists_new):
    p = Point(coord)
    attrs = {'mi': dist * 0.000621371, 'ele': coord[2]}
    f = geojson.Feature(geometry=p, properties=attrs)
    features.append(f)

fc = geojson.FeatureCollection(features)

with open('../hmpoints.geojson', 'w') as f:
    geojson.dump(fc, f)