Hi all,
I am trying to get some mapbox curves into rhino via the mapbox api and rhino3dm library, however I can’t seem to get rid of the tile boundaries which come along with the data. I’ve tried a few methods such as removing the last and first points or setting a tolerance, but none of these works. I wouldn’t even mind a grasshopper solution if there is one.
Thanks
import requests
import math
import mapbox_vector_tile
from pyproj import *
import mercantile
import rhino3dm as rh
lat = -33.8226
lon = 151.2012
zoom = 16
forma = 'mvt'
def deg2num(lat_deg, lon_deg, zoom):
lat_rad = math.radians(lat_deg)
n = 2.0 ** zoom
xtile = int((lon_deg + 180.0) / 360.0 * n)
ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
return (xtile, ytile)
num_x, num_y = deg2num(lat, lon, 16)
def create_boundary(lat, lon, distance):
R = 6371.0
lat_r = math.radians(lat)
ns_dist = distance / (R * 1000)
ew_dist = distance / (R * 1000 * math.cos(lat_r))
max_lat = lat + ns_dist
min_lat = lat - ns_dist
max_lon = lon + ew_dist
min_lon = lon - ew_dist
return min_lon, max_lon, min_lat, max_lat
def create_layer(model, name, color):
layer = rh.Layer()
layer.Name = name
layer.Color = color
return model.Layers.Add(layer)
xmin_LL, xmax_LL, ymin_LL, ymax_LL = create_boundary(lat, lon, 20000)
model = rh.File3dm()
model.Settings.ModelUnitSystem = rh.UnitSystem.Meters
topography_layerIndex = create_layer(model, "Topography", (191, 191, 191, 255))
transformer2 = Transformer.from_crs("EPSG:4326", "EPSG:32756", always_xy=True)
tiles = list(mercantile.tiles(xmin_LL, ymin_LL, xmax_LL, ymax_LL, zooms=16))
for tile in tiles:
tile1 = mercantile.Tile(tile.x, tile.y, 16)
bbox = mercantile.bounds(tile1)
lon1, lat1, lon2, lat2 = bbox
terrain_url = f'https://api.mapbox.com/v4/mapbox.mapbox-terrain-v2/{tile1.z}/{tile1.x}/{tile1.y}.{forma}?access_token={mapbox_access_token}'
terrain_response = requests.get(terrain_url)
terrain_data = terrain_response.content
terrains = mapbox_vector_tile.decode(terrain_data)
terrain_layer = terrains.get('contour')
if terrain_layer:
for feature in terrain_layer['features']:
geometry_type = feature['geometry']['type']
height = feature['properties']['ele']
if geometry_type == 'Polygon':
geometry = feature['geometry']['coordinates']
for ring in geometry:
points = []
for coord in ring:
x_val, y_val = coord[0], coord[1]
x_prop = (x_val / 4096)
y_prop = (y_val / 4096)
lon_delta = lon2 - lon1
lat_delta = lat2 - lat1
lon_mapped = lon1 + (x_prop * lon_delta)
lat_mapped = lat1 + (y_prop * lat_delta)
lon_mapped, lat_mapped = transformer2.transform(
lon_mapped, lat_mapped)
point = rh.Point3d(lon_mapped, lat_mapped, 0)
points.append(point)
polyline = rh.Polyline(points)
curve = polyline.ToNurbsCurve()
att = rh.ObjectAttributes()
att.LayerIndex = topography_layerIndex
att.SetUserString("Elevation Height", str(height))
model.Objects.AddCurve(curve, att)
elif geometry_type == 'MultiPolygon':
geometry = feature['geometry']['coordinates']
for polygon in geometry:
for ring in polygon:
points = []
for coord in ring:
x_val, y_val = coord[0], coord[1]
x_prop = (x_val / 4096)
y_prop = (y_val / 4096)
lon_delta = lon2 - lon1
lat_delta = lat2 - lat1
lon_mapped = lon1 + (x_prop * lon_delta)
lat_mapped = lat1 + (y_prop * lat_delta)
lon_mapped, lat_mapped = transformer2.transform(
lon_mapped, lat_mapped)
point = rh.Point3d(lon_mapped, lat_mapped, 0)
points.append(point)
polyline = rh.Polyline(points)
curve = polyline.ToNurbsCurve()
att = rh.ObjectAttributes()
att.LayerIndex = topography_layerIndex
att.SetUserString("Elevation Height", str(height))
model.Objects.AddCurve(curve, att)
model.Write("C:/Users/RivinduB/Desktop/test.3dm")
test.3dm (375.5 KB)
