diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..bf750d2 --- /dev/null +++ b/.flake8 @@ -0,0 +1,2 @@ +[flake8] +max-line-length = 99 \ No newline at end of file diff --git a/.gitignore b/.gitignore index 0ba19fc..87cb30c 100644 --- a/.gitignore +++ b/.gitignore @@ -161,3 +161,9 @@ cython_debug/ cache +*.stl +junk +proto + + + diff --git a/README.md b/README.md index d942f21..19bfa1f 100644 --- a/README.md +++ b/README.md @@ -1,22 +1,22 @@ -# landscape2stl: High resolution terrain models for 3D printing +# landscape2stl: Modular, high resolution terrain models for 3D printing [Source](https://github.com/gecrooks/landscape2stl) -Gavin E. Crooks (2023) +Gavin E. Crooks (2024) -![Yosemite](images/yosemite_r2d2.jpeg) -Yosemite valley and environs, scale about 1" to 1 mile, overall size 24"x21". R2D2 for scale. +![Yosemite](images/yosemite.jpeg) +Yosemite valley and environs, scale about 1" to 1 mile, overall size 42"x43" (a little over 1 square meter). ![Magnets](images/landscape_magnets.jpeg) -Yosemite was printed in 9 segments, held together by magnets installed along the base. +Yosemite was printed in 30 segments, held together by magnets installed along the base. The python script landscape2stl.py generates large scale terrain models as STL files suitable for 3D printing. I have -incorporated a number of innovations, namely a high resolution data source, novel map projection, and base alignment magnets. +incorporated a number of innovations, notably a high resolution data source and base alignment magnets. Install dependencies, @@ -35,99 +35,105 @@ dataset resolution isn't fine enough the model starts looking like Minecraft.) The downside is that this dataset is only available for the contiguous 48 states, Hawaii, and some parts of Alaska. -## There is no projection - -When you print a 2D map of a portion of the world you have to choose a map projection, since the world is spherical and the map is flat. There are many possible projections, and no projection is ideal since they must all inevitable distort the spherical reality to the flat plane. - -While I was pondering the problem of which projection would be best for terrain models, a flash of inspiration hit me: I don't need a projection at all! With a 3D printer I can print in 3D! Each terrain segment created by landscape2stl.py is a gently curved wedge of the Earth. The sides slope inwards slightly, and the east-west lines of longitude are curved. At large scales this effect is subtle, but not insignificant even at a scale of 1 miles: 1 inch. - ## Base alignment magnets Since print beds are limited in size large landscapes have to be printed as a collection of segments. However, getting neighboring segments to line up satisfactorily can be tricky. Inspired by the gridfinity project, I solved the alignment problem by adding holes along the base into which you can press-fit 6mm X 2mm -magnets. (I'm using DIYMAG brand from Amazon, which cost about 3¢ each.) Neighboring segments snap together with a satisfying click, with nigh on perfect alignment, and can be broken apart again for storage or transport. +magnets. (I'm using DIYMAG brand from Amazon, which cost about 5¢ each.) Neighboring segments snap together with a satisfying click, with nigh on perfect alignment, and can be broken apart again for storage or transport. + +I've also added smaller holes along the sides for alignment pins (Use metal pins, or short lengths of filament). ## CLI: landscape2stl.py ``` -> python -m landscape2stl.py --help +> python -m landscape2stl --help -usage: landscape2stl.py [-h] - [--preset REGION] - [--scale SCALE] [--exaggeration EXAGGERATION] [--magnets MAGNETS] - [N W S E ...] +usage: landscape2stl.py [-h] + [--preset {half_dome,west_of_half_dome,whitney,grand_canyon,shasta,shasta_west,joshua_tree,owens_valley,denali}] [--quad QUAD] + [--state STATE] + [--scale SCALE] + [--exaggeration EXAGGERATION] + [--resolution {10,30}] + [--magnets MAGNETS] + [--name NAME] [-v] + [S W N E ...] -Create landscape STLs +Create quadrangle landscape STLs positional arguments: - N W S E Latitude/longitude coordinates for slab (Order North edge, west edge, south edge, east edge) + S W N E Latitude/longitude coordinates for quadrangle (Order south edge, west edge, north edge, east edge) options: -h, --help show this help message and exit - --preset REGION + --preset {half_dome,west_of_half_dome,whitney,grand_canyon,shasta,shasta_west,joshua_tree,owens_valley,denali} + --quad QUAD + --state STATE --scale SCALE Map scale --exaggeration EXAGGERATION Vertical exaggeration + --resolution {10,30} DEM resolution --magnets MAGNETS Magnet spacing (in degrees) - --name FILENAME - + --name NAME Filename for model + -v, --verbose ``` -### N W S E +### S W N E + +Latitude/longitude coordinates for slab (Order south edge, west edge, north edge, east edge) + -Latitude/longitude coordinates for slab (Order North edge, west edge, south edge, east edge) ### Preset A collection of pre-selected regions that override explicit selection of coordinates. These regions were used for development and illustration. Take a look at the code to see what the options are. +### Quad +Alternatively give the name of a USGS 7.5 minute quadrangle map. This will print at 1:62_500 scale. + +One resource to find the name of quadrangle maps is https://livingatlas.arcgis.com/topomapexplorer/ +(Select 1:24_000 scale for 7.5 minute maps.) + ### Scale Common scales for US maps include: -* 1:24_000, 1" = 2000', about 2.5" to 1 mile -* 1:62_500, about 1" to 1 mile -* 1:125_000, about 1" to 2 miles -* 1:1_000_000, about 1" to 16 miles + +* 1:1_000_000 Approx 1 inch: 16 miles 2° x 2° 144 x 112 miles +* 1:250_000 Approx 1 inch: 4 miles 1/2° x 1/2° +* 1:62_500 Approx 1 inch: 1 mile 1/8° x 1/8° 9 x 7 miles +* 1:24_000 1" = 2000', about 2.5" to 1 mile +* 1:15_625 Approx 4 inches: 1 mile 1/32° x 1/32° 2 1/4 x 1 3/4 miles + +Also listed are the approximate scale in inches per mile, the latitude and longitude extend of the quadrangle that will fit on a 256mm x 256mm print bed, and the approximate scaled size in miles. ### Exaggeration -Low scale models require vertical exaggeration else the landscape looks flat and uninteresting. Exaggeration isn't needed at scales of 1:125_000 or higher. Opinions vary as to appropriate exaggeration, but at 1:1_000_000 an exaggeration of 2:1 looks OK. (Just don't go the way of NASA when showing pictures of Maat Mons) +Low scale models require vertical exaggeration else the landscape looks flat and uninteresting. Exaggeration isn't needed at scales of 1:125_000 or higher. Opinions vary as to appropriate exaggeration, but at 1:1_000_000 an exaggeration of 3:1 looks OK. (Just don't go the way of NASA when showing pictures of Maat Mons) If not set explicitly, the exaggeration is set using a heuristic based on the scale. ### Magnet spacing -The spacing between magnets in degrees. If you use a standard scale a reasonable magnet spacing will be chosen for you. Set to 0.0 -to disable magnets. +The spacing between magnets in degrees. If you use a standard scale a reasonable magnet spacing will be chosen for you. ## Slicing and Printing -* Printer: Ender 3 pro +* Printer: Bambu Lab P1S * Nozzle: 0.4mm -* Layer height: 0.12mm -* Infill: 5% -* Filament: eSUN PLA PRO (PLA+), Cool White (print at 210c, 70c plate temperature) +* Layer height: 0.08mm +* Infill: 8% +* Filament: eSUN PLA PRO (PLA+), Cool White (Sometimes described as Cold White) +* Brim (for bed adhesion) The filament I used is a bright, opaque white that highlights details by throwing valleys and folds into shadow. -With a layer height of 0.12mm contour lines are approximately 25' apart at a scale of 1" : 1 mile. - -Bed adhesion is important since there is a tendency for larger prints to warp upwards at the corners. -Upper layers may self-correct, but the base and magnets will be out of alignment. After some frustration here are -the steps I took to ensure success. +With a layer height of 0.08mm contour lines are approximately 15' apart at a scale of 1" : 1 mile. -1) Carefully level the bed. -2) Print with a brim (easily removable with a deburring tool). -3) Increase plate temperature (I'm currently using 70c) -4) Use a glass build plate and spray adhesive (3DLAC Printer Adhesive). -5) Monitor the first layers, and abort if there is a failure to adhere. +After printing, use a deburring tool to clean up bottom edges, lightly sand the sides, and press-fit the magnets. Make sure the magnet orientations are consistent. I use the same orientation on North and East edges, and flip the orientation for the South and West edges. Make sure to use the same magnet orientations for each additional slab. -After printing, use a deburring tool to clean up bottom edges, lightly sand the sides, and press-fit the magnets. Make sure the magnet orientations are consistent. I use the same orientation on North and East edges, and flip the orientation for the south and west edges. Make sure to use the same magnet orientations for each additional slab. - -To fit the magnets, take a stack of magnets, place over the magnet hole (in the correct orientation), and give the top of the stack a sharp tap with a mallet. If the magnet won't stay in its socket, add a small spot of superglue. +To fit the magnets, take a stack of magnets, place over the magnet hole (in the correct orientation), and give the top of the stack a few lights tap with a mallet. If the magnet won't stay in its socket, add a small spot of superglue. ## Known Issues @@ -136,34 +142,23 @@ To fit the magnets, take a stack of magnets, place over the magnet hole (in the This is an experimental prototype package of beta code for my own amusement. There is no support. Use at your own risk. Caveat emptor. +### Non-manifold models +The generated models have bugs. The Bambu Labs slicer will complain about "Non-manifold edges". Often these errors can be ignored, but sometimes the slicer gets upset, and the model needs to be repaired (Unfortunately Bambu only provides the repair functionality on Windows, not on Mac OS) The Cura slicer is more forgiving. + +Some of these errors are being generated by hacks in my code that are there to get adequate performance out of the Computational Solid Geometry package in ezdxf, and some seem to be being generated by ezdxf itself. ### Limited data coverage The 3DEP 1/3 arc second data set only covers the USA (And not even all of Alaska). Ideally we would fall back to other lower resolution datasets for the rest of the world. +### Projection -### Model arcs - -Large enough models arc upwards due to the curvature of Earth's surface. The magnets are not strong enough to hold the arc themselves, so it may be necessary to add shims to support the middle segments of the arch. - - -### Flat base edges - -The method for adding the magnet holes is an ugly hack. As a simplification, the surface that the magnets are embedded into is flat. That's correct for the North-South edges, but the East-West edges of the terrain actually curve, creating an odd break between the base and the heights. This can be seem in the Denali preset (Since Denali is far north the curvature of lines of longitude is more pronounced). - -In practice, this mismatch between base and heights doesn't seem to be a problem, but the current solution is inelegant and should be rewritten using proper computational geometry. - +I'm using a Lambert conformal conic projection with standard parallels of 33 and 45 degrees. This projection works well for the contiguous United States, but other projections would be better for Hawaii and Alaska, and for other regions of the world. ### Shorelines -Shorelines are generally not visible in the models. - -(A possible fix would be to reduce the height of the ocean by some number of meters. Unfortunately the 3DEP dataset does not make it clear where land ends and the sea begins. You might think everything above sea level is land, but that produces bad looking coast lines. Setting sea level to 1m oddly produces better results, but not ideal. And also some areas of California, notable Death Valley, are below sea level.) - - -### Circular Contours +We drop ocean areas by a small amount to make shorelines more visible in the models. Unfortunately the 3DEP dataset does not make it clear where land ends and the sea begins. You might think that sea is at sea level, and everything above sea level is land, but that produces bad looking coast lines. Setting sea level to 1m oddly produces better results, but not ideal. Also some areas of California, notable Death Valley, are below sea level. And the sea is not always at zero meters (no not obviously good reason). I've added some heuristics to generate shores lines, but they are under-tested. Use at your own risk. -At the center of each terrain segment up is vertical, but at the edges of the segment up is slightly outwards because each segment is curved. This means the contour lines created by the successive layers of the 3D print are not actually horizontal. This is visible if you create a model of a very flat area, such as a patch of ocean. You get circular contour lines centered at the center of the segment. However, this effect isn't noticeable for any interestingly rugged piece of terrain. diff --git a/example_oahu.py b/example_oahu.py new file mode 100644 index 0000000..6dba960 --- /dev/null +++ b/example_oahu.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python + +from landscape2stl import create_stl, STLParameters + + +params = STLParameters( + scale=250_000, + magnet_spacing=0.05, + exaggeration=2, + magnet_holes=True, + bolt_holes=False, + pin_holes=False, + ) + +ee = -158.3 +nn = 21.30 +delta = 0.1 + +# Oahu +for e in range(0, 7): + for n in range(0, 6): + name = f"oahu_250000_{e}{n}" + north = nn + n * delta + east = ee + e * delta + print(name, north, east, north-delta, east+delta) + create_stl(params, (north, east, north-delta, east+delta), name=name, verbose=True) diff --git a/example_yosemite.py b/example_yosemite.py new file mode 100755 index 0000000..b31e6b5 --- /dev/null +++ b/example_yosemite.py @@ -0,0 +1,47 @@ +#!/usr/bin/env python + +from landscape2stl import create_quad_stl + + +# Names of 7.5 minute quadrangles, north to south, blocks east to west. +yosemite_quadrangles = ( + "Kibbie Lake", + "Tiltill Mountain", + "Piute Mountain", + "Matterhorn Peak", + "Dunderberg Peak", + "Lundy", + # "Negit Island", + "Lake Eleanor", + "hetch_hetchy_reservoir", + "ten_lakes", + "falls_ridge", + "tioga_pass", + "mount_dana", + # "lee_vining", + "Ackerson Mountain", + "tamarack_flat", + "yosemite_falls", + "tenaya_lake", + "vogelsang_peak", + "koip_peak", + # "june_lake", + "El Portal", + "el_capitan", + "half_dome", + "merced_peak", + "mount_lyell", + "mount_ritter", + # "mammoth_mountain", + "Buckingham Mountain", + "Wawona", + "Mariposa grove", + "Sing Peak", + "Timber Knob", + "Cattle Mountain", + # "Crystal Crag", +) + + +for name in yosemite_quadrangles: + create_quad_stl(name, "CA", verbose=True) diff --git a/images/yosemite_r2d2.jpeg b/images/yosemite_r2d2.jpeg deleted file mode 100644 index aacc3ee..0000000 Binary files a/images/yosemite_r2d2.jpeg and /dev/null differ diff --git a/landscape2stl.py b/landscape2stl.py index 881136b..5675aa0 100755 --- a/landscape2stl.py +++ b/landscape2stl.py @@ -1,12 +1,12 @@ #!/usr/bin/env python -# Copyright 2023, Gavin E. Crooks and contributors +# Copyright 2023-2024, Gavin E. Crooks # # This source code is licensed under the MIT License # found in the LICENSE file in the root directory of this source tree. # # -# llandscape2stl: High resolution terrain models for 3D printing +# landscape2stl: High resolution terrain models for 3D printing # # See README for more information # @@ -14,92 +14,167 @@ # Gavin E. Crooks 2023 # +# Disclaimer: This is ugly hacked together prototype code. You have been warned. -# Note: +# Note to self: # latitude is north-south # longitude is west-east (long way around) -import py3dep -import numpy as np -import xarray as xr -import argparse, sys - -import numpy as np -from math import pi +import argparse import math - -from stl import mesh # from package numpy-stl import os -from os import path - +import sys +import zipfile +from dataclasses import dataclass +from typing import Optional, Tuple +import ezdxf +import numpy as np +import pandas as pd +import py3dep +import requests +import us +import xarray as xr +from ezdxf.addons.meshex import stl_dumpb +from ezdxf.addons.pycsg import CSG +from ezdxf.render.forms import cylinder_2p +from numpy.typing import ArrayLike from typing_extensions import TypeAlias -from typing import Any -# Many units and coordinate systems. Use various TypeAlias in desperate effort to keep everything straight -LLA: TypeAlias = Any # latitude, longitude, altitude (in meters) coordinates -ENU: TypeAlias = Any # East, North, Up model coordinates (millimeters) +# Many units and coordinate systems. Use TypeAlias's in desperate effort to +# keep everything straight MM: TypeAlias = float # millimeters Meters: TypeAlias = float # meters Degrees: TypeAlias = float +ECEF: TypeAlias = tuple[ + Meters, Meters, Meters +] # Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates +LLA: TypeAlias = tuple[ + Degrees, Degrees, Meters +] # latitude, longitude, altitude (in meters) coordinates +ENU: TypeAlias = tuple[MM, MM, MM] # East, North, Up model coordinates (millimeters) +BBox: TypeAlias = tuple[ + Degrees, Degrees, Degrees, Degrees +] # Geographic bounding box: south, west, north, east + # Various presets featuring varied terrains for testing and development. -presets = { - "half_dome": [(37.76, -119.57, 37.72, -119.52), 24_000], # 37.75°N 119.53 - "west_of_half_dome": [(37.76, -119.62, 37.72, -119.57), 24_000], # 37.75°N 119.53 - "whitney": [(36.60, -118.33, 36.56, -118.28), 24_000], # High point test - "yosemite_west": [(37.80, -119.80, 37.70, -119.65), 62_500], - "yosemite_valley": [(37.80, -119.65, 37.70, -119.50), 62_500], - "clouds_rest": [(37.80, -119.50, 37.70, -119.35), 62_500], - "grand_canyon": [(36.2000, -112.2000, 36.0000, -111.950), 125_000], - "shasta": [(41.5000, -122.3500, 41.3000, -122.0500), 125_000], - "shasta_west": [(41.5000, -122.6500, 41.3000, -122.3500), 125_000], - "joshua_tree": [(34.0000, -116.0000, 33.8, -115.75), 125_000], - "owens_valley": [(38, -120, 36.5, -118), 1_000_000], - "denali": [(63.5000, -152.0000, 62.5000, -150.0000), 1_000_000], +presets: dict[str, tuple[BBox, int]] = { + # "half_dome": ((37.72, -119.57, 37.76, -119.52), 24_000), # 37.75°N 119.53°W + "half_dome": ((37.72, -119.57, 37.77, -119.51), 24_000), # 37.75°N 119.53°W + "west_of_half_dome": ((37.72, -119.62, 37.76, -119.57), 24_000), + "whitney": ((36.56, -118.33, 36.60, -118.28), 24_000), # High point test + "grand_canyon": ((36.0000, -112.2000, 36.2000, -111.9500), 125_000), + "shasta": ((41.3000, -122.3500, 41.5000, -122.0500), 125_000), + "shasta_west": ((41.3000, -122.6500, 41.5000, -122.3500), 125_000), + "joshua_tree": ((33.8, -116.0000, 34.0000, -115.75), 125_000), + "owens_valley": ((36.5, -120, 38, -118), 1_000_000), + "denali": ((62.5000, -152.0000, 63.5000, -150.0000), 1_000_000), } -standard_scales = [ - 24_000, # 1" = 2000', about 2.5" to 1 mile - 62_500, # about 1" to 1 mile - 125_000, # about 1" to 2 miles - 1_000_000, # about 1" to 16 miles -] +# standard_scales = [ +# 24_000, # 1" = 2000', about 2.5" to 1 mile +# 62_500, # about 1" to 1 mile +# 125_000, # about 1" to 2 miles +# 250_000, # about 1" to 4 miles +# 500_000, # about 1" to 8 miles +# 1_000_000, # about 1" to 16 miles +# ] default_cache = "cache" -default_scale = 62_500 -default_steps = 1024 # change to 16 for faster debugging can be helpful -default_resolution: Meters = 10 -meters_to_mm = 1000.0 # TODO: mm_per_meter -default_base: Meters = -100.0 # meters Lowest point in US is -86 m -default_exaggeration = 1.0 -default_magnets: Degrees = 0.05 -default_sea_level: Meters = 1.0 # surprisingly works better than 0 meters - - -default_scale_magnets = { - 24_000: 0.01, - 62_500: 0.025, - 125_000: 0.05, - 1_000_000: 0.25, -} + + +@dataclass +class STLParameters: + """Parameters for the STL terrain models (apart from actual coordinates). + Scale is the important parameter to set, the rest can generally be left at default. + """ + + scale: int = 62_500 + resolution: int = 0 # Auto set in __post_init__ + resolution_choices: tuple[int] = (10, 30) # meters + pitch: MM = 0.40 # Nozzle size + + min_altitude: Meters = -100.0 # Lowest point in US is -86 m + + drop_sea_level: bool = True + sea_level: Meters = 1.7 + sea_level_drop: MM = 0.24 # 3 layers + + exaggeration: float = 0.0 # Auto set in __post_init__ + + base_height: MM = 10.0 + + magnet_holes: bool = True + magnet_spacing: Degrees = 0.0 # Auto set in __post_init__ + magnet_diameter: MM = 6.00 + magnet_padding: MM = 0.025 + magnet_depth: MM = 2.00 + magnet_recess: MM = 0.10 + magnet_sides: int = 24 + + pin_holes: bool = True + pin_length: MM = 9 + pin_diameter: MM = 1.75 + pin_padding: MM = 0.05 * 3 + pin_sides: int = 8 + + bolt_holes: bool = False + bolt_hole_offset: MM = 10 + bolt_hole_diameter: MM = 5.6 + bolt_hole_padding: MM = 0.2 + bolt_hole_depth: MM = 9.1 + bolt_hole_sides: int = 24 + + def __post_init__(self): + if not self.magnet_spacing: + if self.scale == 24_000: + self.magnet_spacing = 1/64 + else: + self.magnet_spacing = self.scale / 2_000_000 + + if not self.resolution: + if self.scale < 250_000: + self.resolution = self.resolution_choices[0] + else: + self.resolution = self.resolution_choices[1] + + if not self.exaggeration: + # Heuristic for vertical exaggeration + # scale exaggeration + # <= 62_500 1.0 + # 125_000 1.5 + # 250_000 2.0 + # 500_000 2.5 + # 1_000_000 3.0 + if self.scale <= 62_500: + self.exaggeration = 1.0 + else: + self.exaggeration = 3 - 0.5 * math.log2(1_000_000/self.scale) + +# end STLParameters def main() -> int: - parser = argparse.ArgumentParser(description="Create landscape STLs") + default_params = STLParameters() + parser = argparse.ArgumentParser(description="Create quadrangle landscape STLs") parser.add_argument( "coordinates", - metavar="N W S E", + metavar="S W N E", type=float, nargs="*", - help="Latitude/longitude coordinates for slab (Order North edge, west edge, south edge, east edge)", + help="Latitude/longitude coordinates for quadrangle (Order south edge, west edge, north edge, east edge)", ) parser.add_argument("--preset", dest="preset", choices=presets.keys()) - parser.add_argument("--scale", dest="scale", type=float, help="Map scale") + parser.add_argument("--quad", dest="quad", type=str) + + parser.add_argument("--state", dest="state", type=str, default="CA") + + parser.add_argument("--scale", dest="scale", type=int, help="Map scale") parser.add_argument( "--exaggeration", @@ -109,12 +184,23 @@ def main() -> int: help="Vertical exaggeration", ) + parser.add_argument( + "--resolution", + dest="resolution", + default=default_params.resolution, + choices=default_params.resolution_choices, + type=int, + help="DEM resolution", + ) + parser.add_argument( "--magnets", dest="magnets", type=float, help="Magnet spacing (in degrees)" ) parser.add_argument("--name", dest="name", type=str, help="Filename for model") + parser.add_argument("-v", "--verbose", action="store_true") + args = vars(parser.parse_args()) name = None @@ -128,250 +214,291 @@ def main() -> int: if args["scale"] is None: args["scale"] = presets[name][1] + if args["name"] is None: + args["name"] = args["preset"] + + if args["quad"] is not None: + name = args["quad"].lower().replace(" ", "_") + coords = quad_coordinates(name, args["state"]) + args["coordinates"] = coords + # args["scale"] = 62_500 # params default to this scale + + if args["name"] is None: + args["name"] = "quad_" + args["state"].lower() + "_" + name + if not args["coordinates"]: parser.print_help() return 0 - if args["magnets"] is None: - args["magnets"] = default_scale_magnets[args["scale"]] + if args["scale"] is None: + args["scale"] = default_params.scale - if args["name"]: - name = args["name"] + params = STLParameters( + scale=args["scale"], + exaggeration=args["exaggeration"], + magnet_spacing=args["magnets"], + resolution=args["resolution"], + # projection=args["projection"], + ) create_stl( - args["coordinates"], - args["scale"], - args["exaggeration"], - magnets=args["magnets"], - name=name, + params, args["coordinates"], filename=args["name"], verbose=args["verbose"] ) return 0 -# todo: rename coords to tile_coords? +def ustopo_current(): + url = "https://prd-tnm.s3.amazonaws.com/StagedProducts/Maps/Metadata/ustopo_current.zip" + zip_file_name = "ustopo_current.zip" + csv_file_name = "ustopo_current.csv" + directory = "cache" + zip_file_path = os.path.join(directory, zip_file_name) -def create_stl( - coords, - scale=default_scale, - exaggeration=default_exaggeration, - steps=default_steps, - resolution=default_resolution, - magnets=default_magnets, - name=None, -): - origin = (coords[0] + coords[2]) / 2, (coords[1] + coords[3]) / 2, -100000 + if not os.path.exists(zip_file_path): + os.makedirs(directory, exist_ok=True) + response = requests.get(url) + with open(zip_file_path, "wb") as f: + f.write(response.content) - # adjust origin so that 4 corners are equal height - north, west, south, east = coords - alt = 0.0 - west_north = lla_to_ecef([north, west, alt]) - east_north = lla_to_ecef([north, east, alt]) - east_south = lla_to_ecef([south, east, alt]) - west_south = lla_to_ecef([south, west, alt]) - new_origin_ecef = [ - (west_north[0] + east_south[0] + west_south[0] + east_north[0]) / 4, - (west_north[1] + east_south[1] + west_south[1] + east_north[1]) / 4, - (west_north[2] + east_south[2] + west_south[2] + east_north[2]) / 4, - ] - new_origin = ecef_to_lla(new_origin_ecef) + with zipfile.ZipFile(zip_file_path, "r") as z: + with z.open(csv_file_name) as csv_file: + df = pd.read_csv(csv_file) - # Honestly, dunno why this works. - origin = origin[0] * 2 - new_origin[0], new_origin[1], 0.0 - crns = corners_to_model(coords, -0, origin, scale, exaggeration) + return df - elevation = download_elevation(coords, steps, resolution) - print("Building terrain...") - surface = elevation_to_surface(elevation, origin, scale, exaggeration) +def quad_coordinates(quad_name, state="CA"): + df = ustopo_current() - # FIXME Ocean - # surface = np.where(surface==0.0, surface-100, surface) - print("Triangulating...") - triangles = triangulate_surface(surface, coords, origin, scale, exaggeration) + quad_name = quad_name.lower().replace("_", " ") - base_triangles = triangulate_base(coords, origin, scale, exaggeration, magnets) + condition = (df["map_name"].str.lower() == quad_name) & df[ + "state_list" + ].str.contains(state) - triangles = np.concatenate([triangles, base_triangles]) + row = df[condition] - min_height = np.min(triangles) + if len(row) == 0: + raise ValueError("Quadrangle " + quad_name + " not found") - if name is None: - filename = "{:.2f}_{:.2f}_{:.2f}_{:.2f}.stl".format(*coords) - else: - filename = name + ".stl" + southbc = row["southbc"].astype(float).iloc[0] + westbc = row["westbc"].astype(float).iloc[0] - data = np.zeros(triangles.shape[0], dtype=mesh.Mesh.dtype) - data["vectors"] = triangles - the_mesh = mesh.Mesh(data.copy()) + return southbc, westbc, southbc + 1 / 8, westbc + 1 / 8 - print(f"Saving {filename}") - the_mesh.save(filename) - print() +def quad_from_coordinates(lat, long): + df = ustopo_current() + condition = ( + (df["southbc"].astype(float) <= lat) + & (lat < df["northbc"].astype(float)) + & (df["westbc"].astype(float) <= long) + & (long < df["eastbc"].astype(float)) + ) -def triangle_normal(A, B, C): - A = np.asarray(A) - B = np.asarray(B) - C = np.asarray(C) + row = df[condition] + print(len(row)) + if len(row) == 0: + return (None, None) + name = row["map_name"].astype(str).iloc[0] + state_name = row["primary_state"].astype(str).iloc[0] + state = us.states.lookup(state_name) - AB = B - A - AC = C - A + return name, state.abbr - normal = np.cross(AB, AC) - normal_normalized = normal / np.linalg.norm(normal) - return normal_normalized +def create_quad_stl(name, state, filename=None, verbose=False): + coords = quad_coordinates(name, state) + if filename is None: + filename = "quad_" + state.lower() + "_" + name.lower().replace(" ", "_") + params = STLParameters( + scale=62_500, + ) -def lla_to_ecef(lat_lon_alt: LLA): - """Convert latitude, longitude, altitude (LLA) coordinates to Earth-Centered, Earth-Fixed (ECEF) Cartesian - coordinates""" - latitude, longitude, altitude = lat_lon_alt + create_stl(params, coords, filename, verbose) - # Constants for WGS84 - a = 6378137.0 # Equatorial radius (meters) - e = 0.08181919084 # Eccentricity - # e = 0 # Approximate earth as sphere - # Convert latitude and longitude to radians - lat_rad = math.radians(latitude) - lon_rad = math.radians(longitude) +def create_stl( + params: STLParameters, + boundary: BBox, + filename: Optional[str] = None, + verbose: bool = False, +) -> None: + # Locate origin + south, west, north, east = boundary + origin = locate_origin(boundary) - # Calculate prime vertical radius of curvature - N = a / math.sqrt(1 - e**2 * math.sin(lat_rad) ** 2) + # Calculate steps... + north_west_enu = lla_to_model((north, west, 0.0), origin, params) + south_east_enu = lla_to_model((south, east, 0.0), origin, params) - # Calculate ECEF coordinates - X = (N + altitude) * math.cos(lat_rad) * math.cos(lon_rad) - Y = (N + altitude) * math.cos(lat_rad) * math.sin(lon_rad) - Z = ((1 - e**2) * N + altitude) * math.sin(lat_rad) + extent_ns = south_east_enu[0] - north_west_enu[0] + extent_we = north_west_enu[1] - south_east_enu[1] - return X, Y, Z + ns_steps = int(round(extent_ns / params.pitch)) + we_steps = int(round(extent_we / params.pitch)) + steps = max(ns_steps, we_steps) + elevation = download_elevation(boundary, steps, params.resolution, verbose) -def ecef_to_lla(ecef): - x, y, z = ecef + if verbose: + print("Building terrain...") - # Constants for WGS84 - a = 6378137.0 # Semi-major axis - e_sq = 0.00669437999014 # Square of eccentricity - # e_sq = 0.0 + surface = elevation_to_surface(elevation, origin, params) - # Calculate longitude - lon = math.atan2(y, x) + # Add a little bit of noise. Hack for smooth seascapes + # (I think this is necessary to make sure we don't have large number of co-planer triangles + # Which seems to upset the CSG module) + # FIXME: Do I still need this hack? + # surface += 0.001 * np.random.uniform(size=surface.shape) - # Iterative calculation of latitude and altitude - p = math.sqrt(x**2 + y**2) - lat = math.atan2(z, p * (1.0 - e_sq)) + if verbose: + print("Triangulating surface...") - while True: - N = a / math.sqrt(1 - e_sq * math.sin(lat) ** 2) - alt = p / math.cos(lat) - N - new_lat = math.atan2(z, p * (1.0 - e_sq * N / (N + alt))) - if abs(new_lat - lat) < 1e-9: - lat = new_lat - break + # Another hack: We build the surface and the base separately. The base + # alone uses Constructive Solid Geometry (CSG). The CSG module in ezdxf is using + # a binary space partitioning (BSP) tree. This shatters triangles into lots of + # sub-triangles. So if we tried using CSG on our landscape surface that the 100s of + # thousands of initial triangles ends up as millions. - lat = new_lat + # Also even with small models ezdxf generates non-manifold meshes that can upset the slicer. + # Errors seem to be reparable or ignorable (mostly). - # Convert from radians to degrees - lon_deg = math.degrees(lon) - lat_deg = math.degrees(lat) + surface_mesh = triangulate_surface(surface, boundary, origin, params) - return lat_deg, lon_deg, alt + if verbose: + print("Triangulating base...") + base_mesh = triangulate_base(boundary, origin, params, steps) -def lla_to_enu(lat, lon, alt, origin_lat, origin_lon, origin_alt): - """Convert latitude, longitude, and altitude (LLA) coordinates to local Cartesian coordinates - (East-North-Up or ENU) relative to a given origin point.""" - # Convert origin to ECEF coordinates - x_target, y_target, z_target = lla_to_ecef([lat, lon, alt]) - x_origin, y_origin, z_origin = lla_to_ecef([origin_lat, origin_lon, origin_alt]) + model = surface_mesh + model.add_mesh(mesh=base_mesh) - # Calculate ECEF vector between origin and target point - dx, dy, dz = x_target - x_origin, y_target - y_origin, z_target - z_origin - # print(dx, dy, dz) + model = model.optimize_vertices() + model.normalize_faces() - # Convert origin latitude and longitude to radians - lat_rad = math.radians(origin_lat) - lon_rad = math.radians(origin_lon) + if verbose: + print("Faces:", len(model.faces)) - # Define the rotation matrix - R = np.array( - [ - [-math.sin(lon_rad), math.cos(lon_rad), 0], - [ - -math.sin(lat_rad) * math.cos(lon_rad), - -math.sin(lat_rad) * math.sin(lon_rad), - math.cos(lat_rad), - ], - [ - math.cos(lat_rad) * math.cos(lon_rad), - math.cos(lat_rad) * math.sin(lon_rad), - math.sin(lat_rad), - ], - ] - ) + if verbose: + print("Creating STL...") - # Multiply the rotation matrix by the ECEF vector - enu = R.dot(np.array([dx, dy, dz])) + if filename is None: + filename = "{:f}_{:f}_{:f}_{:f}.stl".format(*boundary) + else: + filename = filename + ".stl" - return enu[0], enu[1], enu[2] + binary_stl = stl_dumpb(model) + if verbose: + print(f"Saving {filename}") -def lla_to_model( - lat, lon, alt, origin_lat, origin_lon, origin_alt, scale, exaggeration -): - """Convert latitude, longitude, and altitude (LLA) coordinates to model ENU Cartesian coordinates (millimeters)""" - enu = lla_to_enu(lat, lon, alt * exaggeration, origin_lat, origin_lon, origin_alt) - enu = np.asarray(enu) - enu /= scale - # enu *= (1.0, 1.0, exaggeration) - enu *= meters_to_mm + with open(filename, "wb") as binary_file: + binary_file.write(binary_stl) + + +# end create_stl + + +# FIXME: Still needed with current projection? +def locate_origin(boundary: BBox) -> LLA: + # Locate a point in the middle of the landscape tile which defines the + # local up direction. Model distances are defined relative to this point + + south, west, north, east = boundary - return enu + origin = (south + north) / 2, (east + west) / 2, -100000.0 + # adjust origin so that 4 corners are equal height + + alt = -100 + west_north = lla_to_ecef((north, west, alt)) + east_north = lla_to_ecef((north, east, alt)) + east_south = lla_to_ecef((south, east, alt)) + west_south = lla_to_ecef((south, west, alt)) + new_origin_ecef = ( + (west_north[0] + east_south[0] + west_south[0] + east_north[0]) / 4, + (west_north[1] + east_south[1] + west_south[1] + east_north[1]) / 4, + (west_north[2] + east_south[2] + west_south[2] + east_north[2]) / 4, + ) + new_origin = ecef_to_lla(new_origin_ecef) + + # Honestly, dunno why this works. + origin = (origin[0] * 2 - new_origin[0], new_origin[1], 0.0) + + return origin -def download_elevation(coords, steps=default_steps, resolution=default_resolution): - north, west, south, east = coords + +def download_elevation( + boundary: BBox, + steps: int, + resolution: int, + verbose: bool = False, +) -> xr.Dataset | xr.DataArray: + elevation: xr.Dataset | xr.DataArray + south, west, north, east = boundary xcoords = np.linspace(west, east, steps) ycoords = np.linspace(south, north, steps) - filename = "{:.2f}_{:.2f}_{:.2f}_{:.2f}.nc".format(*coords) + filename = "{}_{:.2f}_{:.2f}_{:.2f}_{:.2f}.nc".format(steps, *boundary) if not os.path.exists(default_cache): os.mkdir(default_cache) - fname = os.path.join(default_cache, filename) if not os.path.exists(fname): - print("Downloading elevation data... ", end="", flush=True) + if verbose: + print("Downloading elevation data... ", end="", flush=True) elevation = py3dep.elevation_bygrid( - xcoords, ycoords, crs="EPSG:4326", resolution=10 + tuple(xcoords), tuple(ycoords), crs="EPSG:4326", resolution=resolution ) # or 30 + elevation.to_netcdf(fname) - print("Done", flush=True) + if verbose: + print("Done", flush=True) - print("Loading elevation data from cache...", end="", flush=True) + if verbose: + print("Loading elevation data from cache...", end="", flush=True) elevation = xr.open_dataset(fname) - print("", flush=True) + if verbose: + print("", flush=True) return elevation -def elevation_to_surface(elevation, origin, scale, exaggeration): +def elevation_to_surface( + elevation: xr.Dataset | xr.DataArray, origin: LLA, params: STLParameters +) -> np.ndarray: ycoords = np.asarray(elevation.coords["y"]) xcoords = np.asarray(elevation.coords["x"]) steps = len(ycoords) - elevation = np.asarray(elevation.to_array()).reshape((steps, steps)).T - elevation = np.nan_to_num(elevation, nan=0.0) + elevation_array = np.asarray(elevation.to_array()).reshape((steps, steps)).T - # Uncomment this next line to drop the elevation of ocean - # elevation[elevation <= default_sea_level] = -100 + # Missing date will be nan + elevation_array = np.nan_to_num(elevation_array, nan=0.0) + + if params.drop_sea_level: + from scipy import signal + + dropped_sea_level = ( + -(params.scale * params.sea_level_drop / 1000) / params.exaggeration + ) + sea = np.abs(elevation_array) <= params.sea_level + kernel = np.ones((3, 3), dtype=np.int8) + kernel[1, 1] = 0 + N = signal.convolve(sea, kernel, mode="same") + T = signal.convolve(np.ones_like(elevation_array), kernel, mode="same") + elevation_array = np.where( + np.abs(elevation_array) <= params.sea_level, + dropped_sea_level * (N / T), + elevation_array, + ) surface = np.zeros(shape=(steps, steps, 3)) @@ -379,89 +506,42 @@ def elevation_to_surface(elevation, origin, scale, exaggeration): for y in range(steps): lat = ycoords[y] lon = xcoords[x] - alt = elevation[x, y] - surface[x, y] = lla_to_model(lat, lon, alt, *origin, scale, exaggeration) + alt = elevation_array[x, y] + surface[x, y] = lla_to_model((lat, lon, alt), origin, params) return surface -def normalize(v): - return v / np.linalg.norm(v) - - -def angle_between(v, w): - return np.arccos(np.dot(v, w) / (np.linalg.norm(v) * np.linalg.norm(w))) - - -def rotation_matrix(axis, theta): - kx, ky, kz = axis - c = np.cos(theta) - s = np.sin(theta) - C = 1 - c - - return np.array( - [ - [kx * kx * C + c, kx * ky * C - kz * s, kx * kz * C + ky * s], - [kx * ky * C + kz * s, ky * ky * C + c, ky * kz * C - kx * s], - [kx * kz * C - ky * s, ky * kz * C + kx * s, kz * kz * C + c], - ] - ) - - -def rotate_vector(v, w): - v_normalized = normalize(v) - w_normalized = normalize(w) - - axis = np.cross(v_normalized, w_normalized) - axis_normalized = normalize(axis) - - theta = angle_between(v_normalized, w_normalized) - - R = rotation_matrix(axis_normalized, theta) - - return R - - -def corners_to_model(coords, alt, origin, scale, exaggeration): - north, west, south, east = coords - west_north: ENU = lla_to_model(north, west, alt, *origin, scale, exaggeration) - east_north: ENU = lla_to_model(north, east, alt, *origin, scale, exaggeration) - east_south: ENU = lla_to_model(south, east, alt, *origin, scale, exaggeration) - west_south: ENU = lla_to_model(south, west, alt, *origin, scale, exaggeration) - - return west_north, west_south, east_south, east_north - - -def triangulate_surface(surface, coords, origin, scale, exaggeration): - triangles = [] +def triangulate_surface( + surface: np.ndarray, + boundary: BBox, + origin: LLA, + params: STLParameters, +) -> ezdxf.render.MeshBuilder: + model = ezdxf.render.MeshBuilder() steps = surface.shape[0] - north, west, south, east = coords + south, west, north, east = boundary # Top surface for x in range(steps - 1): for y in range(steps - 1): if ((x + y) % 2) == 0: - triangles.append( + model.add_face( [surface[x, y], surface[x + 1, y], surface[x + 1, y + 1]] ) - triangles.append( + model.add_face( [surface[x, y], surface[x + 1, y + 1], surface[x, y + 1]] ) else: - triangles.append([surface[x, y + 1], surface[x, y], surface[x + 1, y]]) - triangles.append( + model.add_face([surface[x, y + 1], surface[x, y], surface[x + 1, y]]) + model.add_face( [surface[x + 1, y], surface[x + 1, y + 1], surface[x, y + 1]] ) - corners = corners_to_model(coords, default_base, origin, scale, exaggeration) - west_north_bot, west_south_bot, east_south_bot, east_north_bot = corners - bot_west = west_north_bot[0] - bot_north = west_north_bot[1] - bot_alt = west_north_bot[2] - bot_east = east_south_bot[0] - bot_south = east_south_bot[1] + bot_corners = corners_to_model(boundary, params.min_altitude, origin, params) + west_north_bot, west_south_bot, east_south_bot, east_north_bot = bot_corners - bot_height = bot_alt + bot_height = west_north_bot[2] # south xcoords = np.linspace(west_south_bot[0], east_south_bot[0], steps) @@ -469,57 +549,66 @@ def triangulate_surface(surface, coords, origin, scale, exaggeration): for x in range(steps - 1): y = 0 # Sloped edge - tri = ( - surface[x, y], - (xcoords[x], ycoords[x], bot_height), - (xcoords[x + 1], ycoords[x + 1], bot_height), - ) - # curved edge # tri = ( # surface[x, y], - # find_point_on_line(surface[x, y], [surface[x, y][0], surface[x, y][1], default_base], bot_height), - # find_point_on_line(surface[x+1, y], [surface[x+1, y][0], surface[x+1, y][1], default_base], bot_height), - # # (xcoords[x], ycoords[x], bot_height), - # # (xcoords[x + 1], ycoords[x + 1], bot_height), + # (xcoords[x], ycoords[x], bot_height), + # (xcoords[x + 1], ycoords[x + 1], bot_height), # ) - triangles.append(tri) - tri = ( - surface[x, y], - (xcoords[x + 1], ycoords[x + 1], bot_height), - surface[x + 1, y], - ) - # Curved edge + # model.add_face(tri) # tri = ( # surface[x, y], - # find_point_on_line(surface[x+1, y], [surface[x+1, y][0], surface[x+1, y][1], default_base], bot_height), - # # (xcoords[x + 1], ycoords[x + 1], bot_height), + # (xcoords[x + 1], ycoords[x + 1], bot_height), # surface[x + 1, y], # ) - triangles.append(tri) + # model.add_face(tri) - # mag_coords = find_point_on_line(coords0, coords1, mag_height) + tri = ( + surface[x, y], + (surface[x, y][0], surface[x, y][1], bot_height), + (surface[x + 1, y][0], surface[x + 1, y][1], bot_height), + ) + + model.add_face(tri) + tri = ( + surface[x, y], + (surface[x + 1, y][0], surface[x + 1, y][1], bot_height), + surface[x + 1, y], + ) + model.add_face(tri) # north - # print(west_north_bot, east_north_bot) xcoords = np.linspace(west_north_bot[0], east_north_bot[0], steps) ycoords = np.linspace(west_north_bot[1], east_north_bot[1], steps) for x in range(steps - 1): - # print(surface[x, y]) y = -1 + # tri = ( + # surface[x, y], + # (xcoords[x], ycoords[x], bot_height), + # (xcoords[x + 1], ycoords[x + 1], bot_height), + # ) + # model.add_face(tri) + # tri = ( + # surface[x, y], + # (xcoords[x + 1], ycoords[x + 1], bot_height), + # surface[x + 1, y], + # ) + # model.add_face(tri) + tri = ( surface[x, y], - (xcoords[x], ycoords[x], bot_height), - (xcoords[x + 1], ycoords[x + 1], bot_height), + (surface[x, y][0], surface[x, y][1], bot_height), + (surface[x + 1, y][0], surface[x + 1, y][1], bot_height), ) - triangles.append(tri) + + model.add_face(tri) tri = ( surface[x, y], - (xcoords[x + 1], ycoords[x + 1], bot_height), + (surface[x + 1, y][0], surface[x + 1, y][1], bot_height), surface[x + 1, y], ) - triangles.append(tri) + model.add_face(tri) # west xcoords = np.linspace(west_south_bot[0], west_north_bot[0], steps) @@ -531,13 +620,13 @@ def triangulate_surface(surface, coords, origin, scale, exaggeration): (xcoords[s], ycoords[s], bot_height), (xcoords[s + 1], ycoords[s + 1], bot_height), ) - triangles.append(tri) + model.add_face(tri) tri = ( surface[x, s], (xcoords[s + 1], ycoords[s + 1], bot_height), surface[x, s + 1], ) - triangles.append(tri) + model.add_face(tri) # east xcoords = np.linspace(east_south_bot[0], east_north_bot[0], steps) @@ -549,264 +638,488 @@ def triangulate_surface(surface, coords, origin, scale, exaggeration): (xcoords[s], ycoords[s], bot_height), (xcoords[s + 1], ycoords[s + 1], bot_height), ) - triangles.append(tri) + model.add_face(tri) tri = ( surface[x, s], (xcoords[s + 1], ycoords[s + 1], bot_height), surface[x, s + 1], ) - triangles.append(tri) + model.add_face(tri) + + # # Bot surface + # surface[:, :, 2] = bot_height + # for x in range(steps - 1): + # for y in range(steps - 1): + # if ((x + y) % 2) == 0: + # model.add_face( + # [surface[x, y], surface[x + 1, y], surface[x + 1, y + 1]] + # ) + # model.add_face( + # [surface[x, y], surface[x + 1, y + 1], surface[x, y + 1]] + # ) + # else: + # model.add_face([surface[x, y + 1], surface[x, y], surface[x + 1, y]]) + # model.add_face( + # [surface[x + 1, y], surface[x + 1, y + 1], surface[x, y + 1]] + # ) + + return model + + +# FIXME: Needs to be simplified now no longer trying to be overly clever. +def triangulate_base( + boundary: BBox, + origin: LLA, + params: STLParameters, + steps: int, +) -> ezdxf.render.MeshBuilder: + model = ezdxf.render.MeshBuilder() + south, west, north, east = boundary + steps = steps // 8 + + base_alt = params.min_altitude - ( + params.base_height * params.scale / (1000 * params.exaggeration) + ) - # bottom - # triangles.append([west_south_bot, east_south_bot, west_north_bot]) - # triangles.append([east_north_bot, west_north_bot, east_south_bot]) + magnets_alt = (base_alt + params.min_altitude) / 2 - triangles = np.asarray(triangles) - return triangles + top_corners = corners_to_model(boundary, params.min_altitude, origin, params) + west_north_top, west_south_top, east_south_top, east_north_top = top_corners + bot_corners = corners_to_model(boundary, base_alt, origin, params) + west_north_bot, west_south_bot, east_south_bot, east_north_bot = bot_corners -# FIXME: Still not right? -# Magnet bounding rectangle constant altitude rather than z-height + mag_corners = corners_to_model(boundary, magnets_alt, origin, params) + west_north_mag, west_south_mag, east_south_mag, east_north_mag = mag_corners + bot = ( + west_south_bot[2] + east_south_bot[2] + east_north_bot[2] + west_north_bot[2] + ) / 4 -def triangulate_base(coords, origin, scale, exaggeration, magnets): - triangles = [] - north, west, south, east = coords + westing = np.linspace(west, east, steps) + northing = np.linspace(north, south, steps) - default_base_height = 10 # mm - base_alt = default_base - ( - default_base_height * scale / (meters_to_mm * exaggeration) - ) + # South + south_top = [ + lla_to_model((south, w, params.min_altitude), origin, params) for w in westing + ] + south_base = [lla_to_model((south, w, base_alt), origin, params) for w in westing] + south_bot = [ + find_point_on_line(llat, llab, bot) for llat, llab in zip(south_top, south_base) + ] - magnets_alt = (base_alt + default_base) / 2 + for i in range(steps - 1): + model.add_face([south_top[i], south_bot[i], south_top[i + 1]]) + model.add_face([south_top[i + 1], south_bot[i], south_bot[i + 1]]) - top_corners = corners_to_model(coords, default_base, origin, scale, exaggeration) - west_north_top, west_south_top, east_south_top, east_north_top = top_corners + # North + north_top = [ + lla_to_model((north, w, params.min_altitude), origin, params) for w in westing + ] + north_base = [lla_to_model((north, w, base_alt), origin, params) for w in westing] + north_bot = [ + find_point_on_line(llat, llab, bot) for llat, llab in zip(north_top, north_base) + ] - bot_corners = corners_to_model(coords, base_alt, origin, scale, exaggeration) - west_north_bot, west_south_bot, east_south_bot, east_north_bot = bot_corners + for i in range(steps - 1): + model.add_face([north_top[i], north_bot[i], north_top[i + 1]][::-1]) + model.add_face([north_top[i + 1], north_bot[i], north_bot[i + 1]][::-1]) - mag_corners = corners_to_model(coords, magnets_alt, origin, scale, exaggeration) - west_north_mag, west_south_mag, east_south_mag, east_north_mag = mag_corners + # East + east_top = [ + lla_to_model((n, east, params.min_altitude), origin, params) for n in northing + ] + east_base = [lla_to_model((n, east, base_alt), origin, params) for n in northing] + east_bot = [ + find_point_on_line(llat, llab, bot) for llat, llab in zip(east_top, east_base) + ] + + for i in range(steps - 1): + model.add_face([east_top[i], east_bot[i], east_top[i + 1]][::-1]) + model.add_face([east_top[i + 1], east_bot[i], east_bot[i + 1]][::-1]) + + # model.add_face([east_top[0], east_bot[0], east_top[-1]][::-1]) + # model.add_face([east_top[-1], east_bot[0], east_bot[-1]][::-1]) + + # West + west_top = [ + lla_to_model((n, west, params.min_altitude), origin, params) for n in northing + ] + west_base = [lla_to_model((n, west, base_alt), origin, params) for n in northing] + west_bot = [ + find_point_on_line(llat, llab, bot) for llat, llab in zip(west_top, west_base) + ] - # Add base + for i in range(steps - 1): + model.add_face([west_top[i], west_bot[i], west_top[i + 1]]) + model.add_face([west_top[i + 1], west_bot[i], west_bot[i + 1]]) + + # model.add_face([west_top[0], west_bot[0], west_top[-1]]) + # model.add_face([west_top[-1], west_bot[0], west_bot[-1]]) # bot of base - triangles.append([west_south_bot, east_south_bot, west_north_bot]) - triangles.append([east_north_bot, west_north_bot, east_south_bot]) + for i in range(steps - 1): + model.add_face([north_bot[i], south_bot[i], north_bot[i + 1]][::-1]) + model.add_face([north_bot[i + 1], south_bot[i], south_bot[i + 1]][::-1]) # top of base - # triangles.append([west_south_top, east_south_top, west_north_top]) - # triangles.append([east_north_top, west_north_top, east_south_top]) - - bot = west_north_bot[2] - top = west_north_top[2] - mid = west_north_mag[2] - - if magnets == 0.0: - # if True: - # Base sides (No magnets) - triangles.append([east_south_bot, east_south_top, west_south_bot]) - triangles.append([west_south_top, west_south_bot, east_south_top]) - - triangles.append([east_south_top, east_south_bot, east_north_bot]) - triangles.append([east_north_bot, east_north_top, east_south_top]) - - triangles.append([east_north_bot, east_north_top, west_north_bot]) - triangles.append([west_north_top, west_north_bot, east_north_top]) - - triangles.append([west_south_bot, west_south_top, west_north_bot]) - triangles.append([west_north_top, west_north_bot, west_south_top]) - - triangles = np.asarray(triangles) - return triangles - - def trianglate_hole( - mag_coords, mag_normal, scale, exaggeration, mheight=8, mwidth=8 - ): - mt = [] - mheight = 8 - mwidth = 8 - R = rotate_vector((0, 0, 1), mag_normal) - mag_triangles, mag_square = magnet_hole(height=mheight, width=mwidth) - for tri in mag_triangles: - A, B, C = tri - a = R.dot(A) + mag_coords - b = R.dot(B) + mag_coords - c = R.dot(C) + mag_coords - - mt.append([a, b, c]) - - N, W, S, E = mag_square - n = R.dot(N) + mag_coords - w = R.dot(W) + mag_coords - s = R.dot(S) + mag_coords - e = R.dot(E) + mag_coords - - return mt, (n, w, s, e) - - mag_longs = np.arange(west + magnets / 2, east, magnets) - mag_lats = np.arange(south + magnets / 2, north, magnets) + # for i in range(steps - 1): + # model.add_face([north_top[i], south_top[i], north_top[i + 1]]) + # model.add_face([north_top[i + 1], south_top[i], south_top[i + 1]]) + + def make_hole(sides, depth, radius, center, axis): + w = axis + v = (0, 0, 1) + + v_normalized = normalize(v) + w_normalized = normalize(w) + + rot_axis = np.cross(v_normalized, w_normalized) + rot_axis = normalize(rot_axis) + + theta = angle_between(v_normalized, w_normalized) + hole = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + hole.rotate_axis(rot_axis, theta) + hole.translate(*center) + + return hole + + model_csg = CSG(model) + + magnets = params.magnet_spacing long_steps = 1 + round((east - west) / magnets) * 2 - mag2_longs = np.linspace(west, east, long_steps) lat_steps = 1 + round((north - south) / magnets) * 2 - mag2_lats = np.linspace(south, north, lat_steps) - # print(mag2_lats.shape) + longs = np.linspace(west, east, long_steps) + lats = np.linspace(south, north, lat_steps) + + south_normal = triangle_normal(east_south_bot, east_south_top, west_south_bot) + north_normal = triangle_normal(east_north_bot, west_north_bot, east_north_top) + west_normal = triangle_normal(west_south_bot, west_south_top, west_north_top) + east_normal = triangle_normal(east_south_top, east_south_bot, east_north_top) + + if params.magnet_holes: + mag_radius = (params.magnet_diameter) / 2 + params.magnet_padding + mag_depth = params.magnet_depth + params.magnet_recess + mag_sides = params.magnet_sides + + # south + for i in range(1, long_steps, 2): + mag_lla = (south, longs[i], magnets_alt) + mag_enu = lla_to_model(mag_lla, origin, params) + hole = make_hole(mag_sides, mag_depth, mag_radius, mag_enu, south_normal) + model_csg = model_csg - CSG(hole) + + # north + for i in range(1, long_steps, 2): + mag_lla = (north, longs[i], magnets_alt) + mag_enu = lla_to_model(mag_lla, origin, params) + hole = make_hole(mag_sides, mag_depth, mag_radius, mag_enu, north_normal) + model_csg = model_csg - CSG(hole) + + # west + for i in range(1, lat_steps, 2): + mag_lla = (lats[i], west, magnets_alt) + mag_enu = lla_to_model(mag_lla, origin, params) + hole = make_hole(mag_sides, mag_depth, mag_radius, mag_enu, west_normal) + model_csg = model_csg - CSG(hole) + + # east + for i in range(1, lat_steps, 2): + mag_lla = (lats[i], east, magnets_alt) + mag_enu = lla_to_model(mag_lla, origin, params) + hole = make_hole(mag_sides, mag_depth, mag_radius, mag_enu, east_normal) + model_csg = model_csg - CSG(hole) + + if params.pin_holes: + pin_length = params.pin_length + pin_radius = (params.pin_diameter / 2) + params.pin_padding + pin_sides = params.pin_sides + + # south + for i in range(2, long_steps - 1, 2): + pin_lla = (south, longs[i], magnets_alt) + pin_enu = lla_to_model(pin_lla, origin, params) + hole = make_hole(pin_sides, pin_length, pin_radius, pin_enu, south_normal) + model_csg = model_csg - CSG(hole) + + # north + for i in range(2, long_steps - 1, 2): + pin_lla = (north, longs[i], magnets_alt) + pin_enu = lla_to_model(pin_lla, origin, params) + hole = make_hole(pin_sides, pin_length, pin_radius, pin_enu, north_normal) + model_csg = model_csg - CSG(hole) + + # west + for i in range(2, lat_steps - 1, 2): + pin_lla = (lats[i], west, magnets_alt) + pin_enu = lla_to_model(pin_lla, origin, params) + hole = make_hole(pin_sides, pin_length, pin_radius, pin_enu, west_normal) + model_csg = model_csg - CSG(hole) + + # east + for i in range(2, lat_steps - 1, 2): + pin_lla = (lats[i], east, magnets_alt) + pin_enu = lla_to_model(pin_lla, origin, params) + hole = make_hole(pin_sides, pin_length, pin_radius, pin_enu, east_normal) + model_csg = model_csg - CSG(hole) + + if params.bolt_holes: + # Corner bolt holes + offset = params.bolt_hole_offset + sides = params.bolt_hole_sides + radius = params.bolt_hole_padding + params.bolt_hole_diameter / 2 + depth = params.bolt_hole_depth + cylinder = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + cylinder.translate(west_north_bot) + cylinder.translate([offset, -offset, 0]) + model_csg = model_csg - CSG(cylinder) + + cylinder = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + cylinder.translate(west_south_bot) + cylinder.translate([offset, offset, 0]) + model_csg = model_csg - CSG(cylinder) + + cylinder = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + cylinder.translate(east_south_bot) + cylinder.translate([-offset, offset, 0]) + model_csg = model_csg - CSG(cylinder) + + cylinder = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + cylinder.translate(east_north_bot) + cylinder.translate([-offset, -offset, 0]) + model_csg = model_csg - CSG(cylinder) + + center = (np.asarray(east_north_bot) + np.asarray(west_south_bot)) / 2.0 + cylinder = cylinder_2p( + count=sides, + base_center=(0, 0, -depth), + top_center=(0, 0, depth), + radius=radius, + ) + cylinder.translate(center) + model_csg = model_csg - CSG(cylinder) - # print(mag2_lats) + model = model_csg.mesh() + return model - mag_alt = (base_alt + default_base) / 2 - mag_corners = corners_to_model(coords, mag_alt, origin, scale, exaggeration) - mag_height = mag_corners[0][2] - # # south +# End triangulate base - mag_normal = triangle_normal(east_south_bot, east_south_top, west_south_bot) - xcoords_top = np.linspace(west_south_top[0], east_south_top[0], long_steps) - ycoords_top = np.linspace(west_south_top[1], east_south_top[1], long_steps) +def triangle_normal(A: ArrayLike, B: ArrayLike, C: ArrayLike) -> np.ndarray: + A = np.asarray(A) + B = np.asarray(B) + C = np.asarray(C) - xcoords_bot = np.linspace(west_south_bot[0], east_south_bot[0], long_steps) - ycoords_bot = np.linspace(west_south_bot[1], east_south_bot[1], long_steps) + AB = B - A + AC = C - A - xcoords_mag = np.linspace(west_south_mag[0], east_south_mag[0], long_steps) - ycoords_mag = np.linspace(west_south_mag[1], east_south_mag[1], long_steps) + normal = np.cross(AB, AC) + normal_normalized = normal / np.linalg.norm(normal) - for i in range(1, long_steps, 2): - mag_coords = [xcoords_mag[i], ycoords_mag[i], mid] - tri, square = trianglate_hole(mag_coords, mag_normal, scale, exaggeration) - triangles.extend(tri) + return normal_normalized - mag_corners_coords = [ - (xcoords_top[i + 1], ycoords_top[i + 1], top), - (xcoords_top[i - 1], ycoords_top[i - 1], top), - (xcoords_bot[i - 1], ycoords_bot[i - 1], bot), - (xcoords_bot[i + 1], ycoords_bot[i + 1], bot), - ] - for c in [-1, 0, 1, 2]: - tri = square[c], square[c + 1], mag_corners_coords[c] - triangles.append(tri) - tri = square[c + 1], mag_corners_coords[c], mag_corners_coords[c + 1] - triangles.append(tri) - # # north +# FIXME: No longer needed, remove +def lla_to_ecef(lat_lon_alt: LLA) -> ECEF: + """ + Convert latitude, longitude, altitude (LLA) coordinates + to Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates. + """ + latitude, longitude, altitude = lat_lon_alt - mag_normal = triangle_normal(east_north_bot, west_north_bot, east_north_top) + # Constants for WGS84 + a = 6378137.0 # Equatorial radius (meters) + e = 0.08181919084 # Eccentricity - xcoords_top = np.linspace(west_north_top[0], east_north_top[0], long_steps) - ycoords_top = np.linspace(west_north_top[1], east_north_top[1], long_steps) + # Convert latitude and longitude to radians + lat_rad = math.radians(latitude) + lon_rad = math.radians(longitude) - xcoords_bot = np.linspace(west_north_bot[0], east_north_bot[0], long_steps) - ycoords_bot = np.linspace(west_north_bot[1], east_north_bot[1], long_steps) + # Calculate prime vertical radius of curvature + N = a / math.sqrt(1 - e**2 * math.sin(lat_rad) ** 2) - xcoords_mag = np.linspace(west_north_mag[0], east_north_mag[0], long_steps) - ycoords_mag = np.linspace(west_north_mag[1], east_north_mag[1], long_steps) + # Calculate ECEF coordinates + X = (N + altitude) * math.cos(lat_rad) * math.cos(lon_rad) + Y = (N + altitude) * math.cos(lat_rad) * math.sin(lon_rad) + Z = ((1 - e**2) * N + altitude) * math.sin(lat_rad) - for i in range(1, long_steps, 2): - mag_coords = [xcoords_mag[i], ycoords_mag[i], mid] - tri, square = trianglate_hole(mag_coords, mag_normal, scale, exaggeration) - triangles.extend(tri) + return X, Y, Z - mag_corners_coords = [ - (xcoords_top[i + 1], ycoords_top[i + 1], top), - (xcoords_top[i - 1], ycoords_top[i - 1], top), - (xcoords_bot[i - 1], ycoords_bot[i - 1], bot), - (xcoords_bot[i + 1], ycoords_bot[i + 1], bot), - ] - mag_corners_coords = mag_corners_coords[::-1] - for c in [-1, 0, 1, 2]: - tri = square[c], square[c + 1], mag_corners_coords[c] - triangles.append(tri) - tri = square[c + 1], mag_corners_coords[c], mag_corners_coords[c + 1] - triangles.append(tri) +# FIXME: No longer needed, remove +def ecef_to_lla(ecef: ECEF) -> LLA: + """ + Convert to Earth-Centered, Earth-Fixed (ECEF) Cartesian coordinates + to latitude, longitude, altitude (LLA) coordinates. + """ - # west - mag_normal = triangle_normal(west_south_bot, west_south_top, west_north_top) + x, y, z = ecef - mag_lla0 = (mag_lats[0], west, mag_alt) - mag_lla1 = (mag_lats[1], west, mag_alt) - mwidth = ( - lla_to_model(*mag_lla1, *origin, scale, exaggeration)[1] - - lla_to_model(*mag_lla0, *origin, scale, exaggeration)[1] - ) + # Constants for WGS84 + a = 6378137.0 # Semi-major axis + e_sq = 0.00669437999014 # Square of eccentricity - for i in range(1, lat_steps, 2): - # print(i, mag2_lats[i]) - lat = mag2_lats[i] + # Calculate longitude + lon = math.atan2(y, x) - mag_lla0 = (lat, west, base_alt) - mag_lla1 = (lat, west, default_base) - coords0 = lla_to_model(*mag_lla0, *origin, scale, exaggeration) - coords1 = lla_to_model(*mag_lla1, *origin, scale, exaggeration) - mag_coords = find_point_on_line(coords0, coords1, mag_height) + # Iterative calculation of latitude and altitude + p = math.sqrt(x**2 + y**2) + lat = math.atan2(z, p * (1.0 - e_sq)) - tri, square = trianglate_hole( - mag_coords, mag_normal, scale, exaggeration, mwidth=9 - ) - triangles.extend(tri) + while True: + N = a / math.sqrt(1 - e_sq * math.sin(lat) ** 2) + alt = p / math.cos(lat) - N + new_lat = math.atan2(z, p * (1.0 - e_sq * N / (N + alt))) + + if abs(new_lat - lat) < 1e-9: + lat = new_lat + break + + lat = new_lat + + # Convert from radians to degrees + lon_deg = math.degrees(lon) + lat_deg = math.degrees(lat) - xcoords_top = np.linspace(west_south_top[0], west_north_top[0], lat_steps) - ycoords_top = np.linspace(west_south_top[1], west_north_top[1], lat_steps) + return lat_deg, lon_deg, alt - xcoords_bot = np.linspace(west_south_bot[0], west_north_bot[0], lat_steps) - ycoords_bot = np.linspace(west_south_bot[1], west_north_bot[1], lat_steps) - mag_corners_coords = [ - (xcoords_top[i + 1], ycoords_top[i + 1], top), - (xcoords_top[i - 1], ycoords_top[i - 1], top), - (xcoords_bot[i - 1], ycoords_bot[i - 1], bot), - (xcoords_bot[i + 1], ycoords_bot[i + 1], bot), +# FIXME: No longer needed, remove? +def lla_to_enu(lat_lon_alt: LLA, origin_lat_lon_alt: LLA) -> ENU: + """ + Convert latitude, longitude, altitude (LLA) coordinates + to local Cartesian coordinates (East-North-Up or ENU) relative + to a given origin point. In millimeters + """ + # Convert origin to ECEF coordinates + x_target, y_target, z_target = lla_to_ecef(lat_lon_alt) + x_origin, y_origin, z_origin = lla_to_ecef(origin_lat_lon_alt) + + # Calculate ECEF vector between origin and target point + dx, dy, dz = x_target - x_origin, y_target - y_origin, z_target - z_origin + + # Convert origin latitude and longitude to radians + lat_rad = math.radians(origin_lat_lon_alt[0]) + lon_rad = math.radians(origin_lat_lon_alt[1]) + + # Define the rotation matrix + R = np.array( + [ + [-math.sin(lon_rad), math.cos(lon_rad), 0], + [ + -math.sin(lat_rad) * math.cos(lon_rad), + -math.sin(lat_rad) * math.sin(lon_rad), + math.cos(lat_rad), + ], + [ + math.cos(lat_rad) * math.cos(lon_rad), + math.cos(lat_rad) * math.sin(lon_rad), + math.sin(lat_rad), + ], ] - mag_corners_coords = mag_corners_coords[::-1] - for c in [-1, 0, 1, 2]: - tri = square[c], square[c + 1], mag_corners_coords[c] - triangles.append(tri) - tri = square[c + 1], mag_corners_coords[c], mag_corners_coords[c + 1] - triangles.append(tri) + ) - # east - mag_normal = triangle_normal(east_south_top, east_south_bot, east_north_top) + # Multiply the rotation matrix by the ECEF vector + enu = R.dot(np.array([dx, dy, dz])) + enu *= 1000 # meters to mm - mag_lla0 = (mag_lats[0], east, mag_alt) - mag_lla1 = (mag_lats[1], east, mag_alt) - mwidth = ( - lla_to_model(*mag_lla1, *origin, scale, exaggeration)[1] - - lla_to_model(*mag_lla0, *origin, scale, exaggeration)[1] + return enu[0], enu[1], enu[2] + + +def lla_to_model( + lat_lon_alt: LLA, origin_lat_lon_alt: LLA, params: STLParameters +) -> ENU: + """ + Convert latitude, longitude, and altitude (LLA) coordinates + to model ENU Cartesian coordinates in millimeters + """ + + # if params.projection == "none": + # # Probably broken at this point. Do not use. + # assert False + # lat, lon, alt = lat_lon_alt + # enu = lla_to_enu((lat, lon, alt * params.exaggeration), origin_lat_lon_alt) + # enu_scaled = np.asarray(enu) + # enu_scaled /= params.scale + + # return (enu_scaled[0], enu_scaled[1], enu_scaled[2]) + + # elif params.projection == "lambert_conformal_conic": + lat, lon, alt = lat_lon_alt + origin_lat, origin_lon, origin_alt = origin_lat_lon_alt + + east, north = lambert_conformal_conic(lat, lon, center_meridian=origin_lon) + center_east, center_north = lambert_conformal_conic( + origin_lat, origin_lon, center_meridian=origin_lon ) - for i in range(1, lat_steps, 2): - lat = mag2_lats[i] + east = east - center_east + north = north - center_north + alt = alt - origin_alt - mag_lla0 = (lat, east, base_alt) - mag_lla1 = (lat, east, default_base) - coords0 = lla_to_model(*mag_lla0, *origin, scale, exaggeration) - coords1 = lla_to_model(*mag_lla1, *origin, scale, exaggeration) - mag_coords = find_point_on_line(coords0, coords1, mag_height) + up = alt * params.exaggeration + enu_scaled = np.asarray([east, north, up]) + enu_scaled /= params.scale + enu_scaled *= 1000 # meters to mm - tri, square = trianglate_hole(mag_coords, mag_normal, scale, exaggeration) - triangles.extend(tri) + return (enu_scaled[0], enu_scaled[1], enu_scaled[2]) - mag_corners = [ - (mag2_lats[i + 1], east, base_alt), - (mag2_lats[i + 1], east, default_base), - (mag2_lats[i - 1], east, default_base), - (mag2_lats[i - 1], east, base_alt), - ] - mag_corners_coords = [ - lla_to_model(*lla, *origin, scale, exaggeration) for lla in mag_corners - ] - mag_corners_coords = mag_corners_coords[::1] - for c in [-1, 0, 1, 2]: - tri = square[c], square[c + 1], mag_corners_coords[c] - triangles.append(tri) - tri = square[c + 1], mag_corners_coords[c], mag_corners_coords[c + 1] - triangles.append(tri) + # else: + # assert False + + +def normalize(v: np.ndarray) -> np.ndarray: + return v / np.linalg.norm(v) + + +def angle_between(v: np.ndarray, w: np.ndarray) -> np.ndarray: + return np.arccos(np.dot(v, w) / (np.linalg.norm(v) * np.linalg.norm(w))) + + +# FIXME: Probably no longer needed, remove? +def corners_to_model( + boundary: BBox, + alt: Meters, + origin: LLA, + params: STLParameters, +) -> tuple[ENU, ENU, ENU, ENU]: + south, west, north, east = boundary + west_north: ENU = lla_to_model((north, west, alt), origin, params) + east_north: ENU = lla_to_model((north, east, alt), origin, params) + east_south: ENU = lla_to_model((south, east, alt), origin, params) + west_south: ENU = lla_to_model((south, west, alt), origin, params) - triangles = np.asarray(triangles) - return triangles + return west_north, west_south, east_south, east_north +# FIXME: Remove +# This routine no longer necessary now using projection. Should be removed. def find_point_on_line(p1, p2, z3): x1, y1, z1 = p1 x2, y2, z2 = p2 @@ -820,64 +1133,66 @@ def find_point_on_line(p1, p2, z3): return x3, y3, z3 -def magnet_hole(height=10, width=10): - def PointsInCircum(r, n=100, z=0.0): - return [ - (math.cos(2 * pi / n * x) * r, math.sin(2 * pi / n * x) * r, z) - for x in range(0, n + 1) - ] - - magnet_diameter = 6.00 - magnet_depth = 1.65 - magnet_padding = 0.25 - - diameter = 6.5 - depth = 1.8 - sides = 8 - - # hedge = 5 # half edge - top_circle = np.asarray(PointsInCircum(diameter / 2, sides, 0.0)) - bot_circle = np.asarray(PointsInCircum(diameter / 2, sides, -depth)) - center_bot = (0.0, 0.0, -depth) - - triangles = [] - - top_square = [ - [height / 2, 0, 0], - [height / 2, width / 2, 0], - [0, width / 2, 0], - [-height / 2, width / 2, 0], - [-height / 2, 0, 0], - [-height / 2, -width / 2, 0], - [0, -width / 2, 0], - [height / 2, -width / 2, 0], - [height / 2, 0, 0], - ] - - circle = top_circle - for x in range(sides): - tri = (circle[x], circle[x + 1], bot_circle[x]) - triangles.append(tri) - tri = (bot_circle[x + 1], bot_circle[x], circle[x + 1]) - triangles.append(tri) - tri = (bot_circle[x], bot_circle[x + 1], center_bot) - triangles.append(tri) - - tri = ( - circle[x], - top_square[x], - circle[x + 1], - ) - triangles.append(tri) - tri = (top_square[x], top_square[x + 1], circle[x + 1]) - triangles.append(tri) - - return np.asarray(triangles), ( - top_square[1], - top_square[3], - top_square[5], - top_square[7], +def lambert_conformal_conic( + lat: float, + lon: float, + standard_parallel1: float = 33.0, + standard_parallel2: float = 45.0, + center_meridian: float = -96.0, +) -> Tuple[float, float]: + """ + Convert latitude and longitude to Lambert Conformal Conic projection coordinates. + + :param lat: Latitude in degrees. + :param lon: Longitude in degrees. + :param standard_parallel1: First standard parallel. + :param standard_parallel2: Second standard parallel. + :param center_meridian: Longitude of the central meridian. + :return: (x, y) coordinates in the Lambert Conformal Conic projection. + """ + + # Convert degrees to radians + lat = math.radians(lat) + lon = math.radians(lon) + standard_parallel1 = math.radians(standard_parallel1) + standard_parallel2 = math.radians(standard_parallel2) + center_meridian = math.radians(center_meridian) + + # Ellipsoid parameters for WGS 84 + a = 6378137 # semi-major axis + f = 1 / 298.257223563 # flattening + e = math.sqrt(f * (2 - f)) # eccentricity + + # Calculate the scale factor at the standard parallels + m1 = math.cos(standard_parallel1) / math.sqrt( + 1 - e**2 * math.sin(standard_parallel1) ** 2 + ) + m2 = math.cos(standard_parallel2) / math.sqrt( + 1 - e**2 * math.sin(standard_parallel2) ** 2 ) + t = math.tan(math.pi / 4 - lat / 2) / ( + (1 - e * math.sin(lat)) / (1 + e * math.sin(lat)) + ) ** (e / 2) + t1 = math.tan(math.pi / 4 - standard_parallel1 / 2) / ( + (1 - e * math.sin(standard_parallel1)) / (1 + e * math.sin(standard_parallel1)) + ) ** (e / 2) + t2 = math.tan(math.pi / 4 - standard_parallel2 / 2) / ( + (1 - e * math.sin(standard_parallel2)) / (1 + e * math.sin(standard_parallel2)) + ) ** (e / 2) + + # Calculate the scale factor n + n = math.log(m1 / m2) / math.log(t1 / t2) + + # Calculate the projection constants F and rho0 + F = m1 / (n * t1**n) + rho = a * F * t**n + rho0 = a * F * t1**n + + # Calculate the projected coordinates + x = rho * math.sin(n * (lon - center_meridian)) + y = rho0 - rho * math.cos(n * (lon - center_meridian)) + + return x, y if __name__ == "__main__": diff --git a/magnets.stl b/magnets.stl deleted file mode 100644 index c45c008..0000000 Binary files a/magnets.stl and /dev/null differ diff --git a/magnets2stl.py b/magnets2stl.py index 4fdeaaf..49c3fe0 100755 --- a/magnets2stl.py +++ b/magnets2stl.py @@ -19,106 +19,92 @@ # Gavin E. Crooks 2023 -import math -from math import pi -import numpy as np -from stl import mesh +import ezdxf +from ezdxf.render.forms import cube, cylinder_2p +from ezdxf.addons.pycsg import CSG +from ezdxf.addons.meshex import stl_dumps -magnet_diameter = 5.95 # magnet diameter in mm (measured with calipers) -magnet_depth = 1.65 # magnet depth in mm (measured with calipers) -magnet_padding = 0.25 +# create new DXF document +doc = ezdxf.new() +msp = doc.modelspace() +from landscape2stl import STLParameters -diameter = magnet_diameter + magnet_padding * 2 -depth = magnet_depth + magnet_padding -sides = 8 # Magnet hole is octagonal prism. Seems to work fine. +model = ezdxf.render.MeshBuilder() -def PointsInCircum(r, n=100, z=0.0): - return [ - (math.cos(2 * pi / n * x) * r, math.sin(2 * pi / n * x) * r, z) - for x in range(0, n + 1) - ] +# create same geometric primitives as MeshTransformer() objects +cube1 = cube() +cube1 = cube1.scale(20, 20, 10) -circle = np.asarray(PointsInCircum(diameter / 2, sides, 0.0)) +model_csg = CSG(cube1) -bot_circle = np.asarray(PointsInCircum(diameter / 2, sides, -depth)) -center_bot = (0.0, 0.0, -depth) +params = STLParameters() +magnet_radius = (params.magnet_diameter)/2 + params.magnet_padding +magnet_depth = params.magnet_depth + params.magnet_recess +magnet_sides = params.magnet_sides -triangles = [] +pin_length = params.pin_length +pin_radius = (params.pin_diameter/2) + params.pin_padding +pin_sides = params.pin_sides -hedge = 5 -top_square = [ - [hedge, 0, 0], - [hedge, hedge, 0], - [0, hedge, 0], - [-hedge, hedge, 0], - [-hedge, 0, 0], - [-hedge, -hedge, 0], - [0, -hedge, 0], - [hedge, -hedge, 0], - [hedge, 0, 0], -] +cylinder = cylinder_2p( + count=magnet_sides, + base_center=(0, -magnet_depth, 0), + top_center=(0, magnet_depth, 0), + radius=magnet_radius, +) +cylinder.translate(-5, 10, 0) +model_csg = model_csg - CSG(cylinder) +cylinder = cylinder_2p( + count=magnet_sides, + base_center=(0, -magnet_depth, 0), + top_center=(0, magnet_depth, 0), + radius=magnet_radius, +) +cylinder.translate(5, 10, 0) +model_csg = model_csg - CSG(cylinder) -for x in range(sides): - tri = (circle[x], circle[x + 1], bot_circle[x]) - triangles.append(tri) - tri = (bot_circle[x + 1], bot_circle[x], circle[x + 1]) - triangles.append(tri) - tri = (bot_circle[x], bot_circle[x + 1], center_bot) - triangles.append(tri) - tri = ( - circle[x], - top_square[x], - circle[x + 1], - ) - triangles.append(tri) - tri = (top_square[x], top_square[x + 1], circle[x + 1]) - triangles.append(tri) +cylinder = cylinder_2p( + count=pin_sides, + base_center=(0, -pin_length, 0), + top_center=(0, pin_length, 0), + radius=1, +) +cylinder.translate(-5,-10, 0) +model_csg = model_csg - CSG(cylinder) +cylinder = cylinder_2p( + count=pin_sides, + base_center=(0, -pin_length, 0), + top_center=(0, pin_length, 0), + radius=1, +) +cylinder.translate(5, -10, 0) +model_csg = model_csg - CSG(cylinder) -tnw = [hedge, hedge, 0] -tne = [hedge, -hedge, 0] -tse = [-hedge, -hedge, 0] -tsw = [-hedge, hedge, 0] -bnw = [hedge, hedge, -hedge] -bne = [hedge, -hedge, -hedge] -bse = [-hedge, -hedge, -hedge] -bsw = [-hedge, hedge, -hedge] -# north -triangles.append([tnw, tne, bne]) -triangles.append([bne, bnw, tnw]) -# triangles.append([tnw, tne, bne]) +offset = params.bolt_hole_offset +sides = params.bolt_hole_sides +radius = params.bolt_hole_padding + params.bolt_hole_diameter/2 +depth = params.bolt_hole_depth +cylinder = cylinder_2p( + count=sides, base_center=(0, 0, -depth), top_center=(0, 0, depth), radius=radius +) +cylinder.translate([0, 0, -5]) +model_csg = model_csg - CSG(cylinder) -# west -triangles.append([tsw, tnw, bsw]) -triangles.append([bnw, bsw, tnw]) -# south -triangles.append([tsw, tse, bsw]) -triangles.append([bse, bsw, tse]) -# east -triangles.append([tne, tse, bse]) -triangles.append([bse, bne, tne]) +s = stl_dumps(model_csg.mesh()) +filename = 'magnets.stl' -# bot -triangles.append([bnw, bne, bse]) -triangles.append([bse, bsw, bnw]) - - -triangles = np.asarray(triangles) - -data = np.zeros(triangles.shape[0], dtype=mesh.Mesh.dtype) -data["vectors"] = triangles -name = "magnets" -the_mesh = mesh.Mesh(data.copy()) -the_mesh.save(name + ".stl") +with open(filename, "w") as f: + f.write(s) diff --git a/requirements.txt b/requirements.txt index 9242b29..565fe48 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,8 @@ numpy scipy -numpy-stl xarray -py3dep \ No newline at end of file +py3dep +ezdxf +pandas +requests +us