From 546d6700d3ed55a275d80bdd7d39ccb165f53bee Mon Sep 17 00:00:00 2001 From: tf17270 Date: Thu, 9 May 2024 16:56:24 +0100 Subject: [PATCH] o3d fix 07 and tidy --- .../07_aperture_farfield_polarisation.py | 2 - .../09_satellite_channel_modelling.py | 377 ------------------ lyceanem/geometry/geometryfunctions.py | 9 - lyceanem/models/frequency_domain.py | 3 +- 4 files changed, 1 insertion(+), 390 deletions(-) delete mode 100644 docs/examples/09_satellite_channel_modelling.py diff --git a/docs/examples/07_aperture_farfield_polarisation.py b/docs/examples/07_aperture_farfield_polarisation.py index c01669b..cfcfb6e 100644 --- a/docs/examples/07_aperture_farfield_polarisation.py +++ b/docs/examples/07_aperture_farfield_polarisation.py @@ -38,7 +38,6 @@ horn_antenna=antenna_structures(structures(solids=[structure]), points(points=[array_points])) -horn_antenna.visualise_antenna() from lyceanem.models.frequency_domain import calculate_farfield @@ -128,7 +127,6 @@ r=R.from_euler('xyz', np.radians(np.asarray([45.0,45.0,0.0]))) horn_antenna.rotate_antenna(r.as_matrix()) -horn_antenna.visualise_antenna() desired_E_axis = np.zeros((1, 3), dtype=np.complex64) diff --git a/docs/examples/09_satellite_channel_modelling.py b/docs/examples/09_satellite_channel_modelling.py deleted file mode 100644 index 14a656c..0000000 --- a/docs/examples/09_satellite_channel_modelling.py +++ /dev/null @@ -1,377 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -Satellite Channel Modelling -====================================================== - -This example uses the frequency domain :func:`lyceanem.models.frequency_domain.calculate_scattering` function to -predict the scattering parameters for the frequency and environment included in the model. This allows the effects of antenna polarisation, - - -""" - -# %% -# Define Map and Ground Station Location and list of Satellites using skyfield API to load starlink satellites. -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - -import folium - -interest_center = [51.458729, -2.602153] # location of Bristol Physics Observatory - -m = folium.Map( - location=[ - interest_center[0], - interest_center[1], - ], - zoom_start=6, - tiles="Stamen Terrain", -) - -import numpy as np - -from skyfield.api import load, wgs84 - -BristolPhysicsObservatory = wgs84.latlon(interest_center[0], interest_center[1], 0) - -receiver = BristolPhysicsObservatory - -# Build the time range `t` over which to plot, plus other values. - -ts = load.timescale() - -time_points = 2 * 60 # track over 10 minutes -# change over time -t = ts.utc(2022, 12, 15, 13, 4, range(0, time_points)) -timebase = t.utc_datetime() - -# %% -# GeoJSON -# ~~~~~~~~ -# The Folium plugin `TimestampedGeoJson`_ will be used to plot our tracks using timestamped -# GeoJSONs. As a result, we want to convert our data into `GeoJSON format`_. Firstly, we create -# our feature collection which we will append our features to. -# -# .. _GeoJSON format: https://geojson.org/ -# .. _TimestampedGeoJson: https://python-visualization.github.io/folium/plugins.html - -geo_features = list() -geo_json = { - "type": "FeatureCollection", - "features": geo_features, -} -# %% -# Plotting Tracks of Satellites visible to Physics Observatory -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# Now we append our tracks to our feature collection list. We define `colour_iter` which will -# allow us to cycle through different colours when mapping our tracks. - -from collections import defaultdict -from itertools import cycle - -colour_iter = iter( - cycle( - [ - "red", - "blue", - "green", - "purple", - "orange", - "darkred", - "#0909FF", - "#F70D1A", - "#FF6700", - "lightgreen", - "#0AFFFF", - "#12AD2B", - "#E2F516", - "#FFFF00", - "#F52887", - ] - ) -) -colour = defaultdict(lambda: next(colour_iter)) - -trail_Size = 7 # trail_Size is the number of timestamps we want track to trail for - - -# %% -# Define Starlinks Tracks -# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -# Convinience function to limit the tracks to a set number of the nearest N satellites. - -def starlink_tracks(pointofinterest, timebase, num_satellites=5): - starlink_url = "https://celestrak.com/NORAD/elements/supplemental/starlink.txt" - starlinks = load.tle_file(starlink_url) - tracks = [] - for sat in starlinks: - pos = sat.at(timebase) - tracks.append(pos) - latcheck = np.any(np.isnan(wgs84.subpoint_of(pos).latitude.degrees)) - loncheck = np.any(np.isnan(wgs84.subpoint_of(pos).longitude.degrees)) - if latcheck | loncheck: - tracks.pop(-1) - print("caught one") - - distance_trim = np.full((len(tracks), len(timebase)), np.inf) - for track in range(len(tracks)): - lat = wgs84.subpoint_of(tracks[track]).latitude - lon = wgs84.subpoint_of(tracks[track]).longitude - relative_measure = (starlinks[track] - pointofinterest).at(timebase) - alt, az, distance1 = relative_measure.altaz() - altitude = alt.degrees - azimuth = az.degrees - relative_range = distance1.km - plot_points = [] - plot_times = [] - az_trim = [] - alt_trim = [] - range_trim = [] - for point in range(len(lat.degrees)): - if altitude[point] >= 0.0: - distance_trim[track, point] = relative_range[point] - - def get_indices_of_k_smallest(arr, k): - idx = np.argpartition(arr.ravel(), k) - return tuple(np.array(np.unravel_index(idx, arr.shape))[:, range(min(k, 0), max(k, 0))]) - - - indexholder = np.zeros((num_satellites, len(timebase)), dtype=int) - for point in range(len(timebase)): - indexholder[:, point] = get_indices_of_k_smallest(distance_trim[:, point], num_satellites)[0] - - return starlinks, indexholder - -starlinks, indexholder = starlink_tracks(BristolPhysicsObservatory, t, num_satellites=20) - -tracks = [] -for sat in starlinks: - pos = sat.at(t) - tracks.append(pos) - latcheck = np.any(np.isnan(wgs84.subpoint_of(pos).latitude.degrees)) - loncheck = np.any(np.isnan(wgs84.subpoint_of(pos).longitude.degrees)) - if latcheck | loncheck: - tracks.pop(-1) - print("caught one") - -for track in range(len(tracks)): - lat = wgs84.subpoint_of(tracks[track]).latitude - lon = wgs84.subpoint_of(tracks[track]).longitude - relative_measure = (starlinks[track] - BristolPhysicsObservatory).at(t) - alt, az, distance = relative_measure.altaz() - altitude = alt.degrees - azimuth = az.degrees - relative_range = distance.km - plot_points = [] - plot_times = [] - az_trim = [] - alt_trim = [] - range_trim = [] - for point in range(len(lat.degrees)): - if altitude[point] >= 0.0 and np.any(np.isin(track, indexholder[:, point])): - plot_points.append((lon.degrees[point], lat.degrees[point])) - plot_times.append(timebase[point].strftime("%Y-%m-%d %H:%M:%S")) - az_trim.append(azimuth[point]) - alt_trim.append(altitude[point]) - range_trim.append(relative_range[point]) - - for time_index, time in enumerate(plot_times): - geo_features.append( - { - "type": "Feature", - "properties": { - "name": starlinks[track].name, - "style": {"color": colour[track], "weight": 6}, - "times": [time] * len(plot_points[: time_index + 1][-trail_Size:]), - }, - "geometry": { - "type": "LineString", - "coordinates": plot_points[: time_index + 1][-trail_Size:], - }, - } - ) - geo_features.append( - { - "type": "Feature", - "properties": { - "icon": "marker", - "iconstyle": { - "iconUrl": f"http://icons.iconarchive.com/icons/google/noto-emoji-travel-places/1024/42597-satellite-icon.png", - "iconSize": [24, 24], - "fillOpacity": 1, - "popupAnchor": [1, -17], - }, - "popup": "Satellite: " + starlinks[track].name + "
" - "Latitude: " - + "%s" % float("%.8g" % plot_points[time_index][0]) - + "
" - "Longitude: " - + "%s" % float("%.8g" % plot_points[time_index][1]) - + "
" # rounding 8 sigfigs - "Azimuth: " - + "%s" % float("%.8g" % az_trim[time_index]) - + "°" - + "
" - "Altitude: " - + "%s" % float("%.8g" % alt_trim[time_index]) - + "°" - + "
" - "Slant Range: " - + "%s" % float("%.8g" % range_trim[time_index]) - + "km", - "name": starlinks[track].name, - "style": {"color": "black", "weight": 2}, - "times": [time], - }, - "geometry": { - "type": "MultiPoint", - "coordinates": [plot_points[time_index]], - }, - } - ) - -# plot receiver -r_track = [] -for time in range(len(t)): - pos = BristolPhysicsObservatory.at(t) - r_track.append(pos) - -for track in range(len(r_track)): - lat = wgs84.subpoint_of(r_track[track]).latitude - lon = wgs84.subpoint_of(r_track[track]).longitude - relative_measure = (receiver - BristolPhysicsObservatory).at(t) - alt, az, distance = relative_measure.altaz() - altitude = alt.degrees - azimuth = az.degrees - relative_range = distance.km - plot_points = [] - plot_times = [] - az_trim = [] - alt_trim = [] - range_trim = [] - for point in range(len(lat.degrees)): - # if altitude[point] >= 0.0: - plot_points.append((lon.degrees[point], lat.degrees[point])) - plot_times.append(timebase[point].strftime("%Y-%m-%d %H:%M:%S")) - az_trim.append(azimuth[point]) - alt_trim.append(altitude[point]) - range_trim.append(relative_range[point]) - - for time_index, time in enumerate(plot_times): - geo_features.append( - { - "type": "Feature", - "properties": { - "name": "Ground Station at Bristol Physics Observatory", - "style": {"color": colour[track], "weight": 6}, - "times": [time] * len(plot_points[: time_index + 1][-trail_Size:]), - }, - "geometry": { - "type": "LineString", - "coordinates": plot_points[: time_index + 1][-trail_Size:], - }, - } - ) - geo_features.append( - { - "type": "Feature", - "properties": { - "icon": "marker", - "iconstyle": { - "iconUrl": f"https://www.bristol.ac.uk/media-library/sites/physics/new-images/501_10013623.jpg", - "iconSize": [64, 32], - "fillOpacity": 1, - "popupAnchor": [1, -17], - }, - "popup": "Ground Station at Bristol Physics Observatory" + "
" - "Latitude: " - + "%s" % float("%.8g" % plot_points[time_index][0]) - + "
" - "Longitude: " - + "%s" % float("%.8g" % plot_points[time_index][1]) - + "
" # rounding 8 sigfigs - "Azimuth: " - + "%s" % float("%.8g" % az_trim[time_index]) - + "°" - + "
" - "Altitude: " - + "%s" % float("%.8g" % alt_trim[time_index]) - + "°" - + "
" - "Slant Range: " - + "%s" % float("%.8g" % range_trim[time_index]) - + "km", - "name": "Ground Station at Bristol Physics Observatory", - "style": {"color": "black", "weight": 2}, - "times": [time], - }, - "geometry": { - "type": "MultiPoint", - "coordinates": [plot_points[time_index]], - }, - } - ) - -from folium.plugins import TimestampedGeoJson, Fullscreen - -Fullscreen().add_to(m) - -import numpy as np -from geopandas import GeoDataFrame -from shapely.geometry import Polygon, MultiPolygon - - -def collec_to_gdf(collec_poly): - """Transform a `matplotlib.contour.QuadContourSet` to a GeoDataFrame""" - polygons, colors = [], [] - for i, polygon in enumerate(collec_poly.collections): - mpoly = [] - for path in polygon.get_paths(): - try: - path.should_simplify = False - poly = path.to_polygons() - # Each polygon should contain an exterior ring + maybe hole(s): - exterior, holes = [], [] - if len(poly) > 0 and len(poly[0]) > 3: - # The first of the list is the exterior ring : - exterior = poly[0] - # Other(s) are hole(s): - if len(poly) > 1: - holes = [h for h in poly[1:] if len(h) > 3] - mpoly.append(Polygon(exterior, holes)) - except: - print('Warning: Geometry error when making polygon #{}' - .format(i)) - if len(mpoly) > 1: - mpoly = MultiPolygon(mpoly) - polygons.append(mpoly) - colors.append(polygon.get_facecolor().tolist()[0]) - elif len(mpoly) == 1: - polygons.append(mpoly[0]) - colors.append(polygon.get_facecolor().tolist()[0]) - return GeoDataFrame( - geometry=polygons, - data={'RGBA': colors}, - crs={'init': 'epsg:4326'}) - -# %% -# The Results -# ~~~~~~~~~~~~~~ - -TimestampedGeoJson( - data=geo_json, - transition_time=200, - auto_play=True, - add_last_point=False, - period="PT1S", - duration="PT0S", -).add_to(m) - -# %% - - -m.show_in_browser() - - - - - diff --git a/lyceanem/geometry/geometryfunctions.py b/lyceanem/geometry/geometryfunctions.py index 7c49122..00fd15c 100644 --- a/lyceanem/geometry/geometryfunctions.py +++ b/lyceanem/geometry/geometryfunctions.py @@ -52,11 +52,7 @@ def mesh_rotate(mesh, rotation, rotation_centre=np.zeros((1, 3), dtype=np.float3 else: raise ValueError("Rotation must be a 3x1 or 3x3 array") if rotation_centre.shape == (3,) or rotation_centre.shape == (3,1): - print("Warning: rotation_centre is a 1x3 array, reshaping to 3x1") rotation_centre = rotation_centre.reshape(1, 3) - if rotation_centre.shape != (1, 3): - print(rotation_centre.shape) - raise ValueError("Rotation centre must be a 3x1 array") rotated_points = r.apply(mesh.points - rotation_centre) + rotation_centre cell_data = mesh.cell_data point_data = mesh.point_data @@ -116,18 +112,13 @@ def mesh_conversion(conversion_object): triangles = RF.convertTriangles(conversion_object) elif isinstance(conversion_object, list): triangles = np.empty((0), dtype=base_types.triangle_t) - # print("Detected List") for item in conversion_object: - # print(type(item)) - # print(isinstance(item,type(o3d.geometry.TriangleMesh()))) if isinstance(item, meshio.Mesh): triangles = np.append(triangles, RF.convertTriangles(item), axis=0) - # print(len(triangles)) elif isinstance(item, base_classes.structures): triangles = np.append( triangles, item.triangles_base_raycaster(), axis=0 ) - # print(len("B ",triangles)) else: print("no structures") diff --git a/lyceanem/models/frequency_domain.py b/lyceanem/models/frequency_domain.py index 920b218..c80177e 100644 --- a/lyceanem/models/frequency_domain.py +++ b/lyceanem/models/frequency_domain.py @@ -511,8 +511,7 @@ def calculate_scattering( if scattering == 0: # only use the aperture point cloud, no scattering required. - scatter_points = o3d.geometry.PointCloud() - + scatter_points = meshio.Mesh(points= np.empty((0, 3)), cells=[]) unified_model = np.append( np.asarray(aperture_coords.points).astype(np.float32), np.asarray(sink_coords.points).astype(np.float32),