Simple survey creation

Kaufmann, 2022-2024.

[1]:

from geometron.survey import TopoPoint, TopoLine
from geometron.geometries import features_to_gdf
from geometron.geometries import transforms as ggt
from geometron.plot import plot_point, plot_line, plot_profile, plot_gdf_survey, plot_basemap
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point, LineString
import shapely

Create some landmarks in a local coordinate system

Coordinates are planar coordinates

Using wkt

[2]:

L001 = TopoPoint('POINT (0. 0.)', id='L001', kind='landmark')
L002 = TopoPoint('POINT (50. 0.)', id='L002', kind='landmark')

Using shapely geometries

[3]:

L003 = TopoPoint(Point(4., 120.), id='L003', kind='landmark')
L004 = TopoPoint(Point(54., 112.), id='L004', kind='landmark')

Plot L001 in a matplotlib axes

[4]:

L001.plot()

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_8_0.png

Add the landmarks in a survey geodataframe and display the survey

[5]:

gdf = features_to_gdf([L001, L002, L003, L004])
gdf

[5]:

	geometry	id	label	kind	class
0	POINT (0.000 0.000)	L001	L001	landmark	TopoPoint
1	POINT (50.000 0.000)	L002	L002	landmark	TopoPoint
2	POINT (4.000 120.000)	L003	L003	landmark	TopoPoint
3	POINT (54.000 112.000)	L004	L004	landmark	TopoPoint

[6]:

ax = plot_gdf_survey(gdf)

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_11_0.png

Create a baseline in local coordinates

Baseline B01 joins L001 and L002

[7]:

B01 = TopoLine(LineString([L001.geometry, L002.geometry]), id='B01', kind='baseline', show_label=False)

[8]:

B01.plot()

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_15_0.png

Baseline B02 joins L003 and L004

[9]:

B02 = TopoLine(LineString([L003.geometry, L004.geometry]), id='B02', kind='baseline', show_label=False)

Append the survey geodataframe

[10]:

gdf = pd.concat([gdf, features_to_gdf([B01, B02])], ignore_index=True)

Display the appended survey

[11]:

plot_gdf_survey(gdf)

[11]:

<Axes: >

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_21_1.png

Add 18 meandering profiles starting and ending on the baselines

The distance between the starting and ending points along the baseline is set to 2 meters

[12]:

number_of_profiles = 18
spacing_on_B01 = 2
spacing_on_B02 = 2

Create the stakes where the profiles start and end

[13]:

stakes_B01 = [TopoPoint(Point(B01.geometry.interpolate(i*spacing_on_B01)), id='', kind='stake') for i in range(1, number_of_profiles+1)]
stakes_B02 = [TopoPoint(Point(B02.geometry.interpolate(i*spacing_on_B02)), id='', kind='stake') for i in range(1, number_of_profiles+1)]

Append the survey geodataframe and plot it

[14]:

gdf = pd.concat([gdf, features_to_gdf(stakes_B01), features_to_gdf(stakes_B02)], ignore_index=True)

[15]:

plot_gdf_survey(gdf)

[15]:

<Axes: >

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_28_1.png

Create the profiles and append the survey geodataframe

[16]:

profiles = []
for i in range(1,number_of_profiles+1,2):
    profiles.append(TopoLine(LineString([stakes_B01[i-1].geometry, stakes_B02[i-1].geometry]), id=f'P{i:02d}', kind='profile'))
for i in range(2,number_of_profiles+1,2):
    profiles.append(TopoLine(LineString([stakes_B02[i-1].geometry, stakes_B01[i-1].geometry]), id=f'P{i:02d}', kind='profile'))

[17]:

gdf = pd.concat([gdf, features_to_gdf(profiles)], ignore_index=True)

[18]:

fig, ax = plt.subplots(figsize=(6,12))
plot_gdf_survey(gdf, ax=ax)

[18]:

<Axes: >

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_32_1.png

Transform from local coordinates to geographic coordinates

Here an affine transform is used. Based on the use case other approaches could be better suited.

Define the control points for the affine transform

[19]:

origin_coords = [L001, L002, L003, L004]

[20]:

landmarks = gpd.read_file('./data/landmarks.shp')
landmarks.set_index('Label', drop=False, inplace=True)
landmarks.crs

[20]:

<Projected CRS: EPSG:3812>
Name: ETRS89 / Belgian Lambert 2008
Axis Info [cartesian]:
- X[east]: Easting (metre)
- Y[north]: Northing (metre)
Area of Use:
- name: Belgium - onshore.
- bounds: (2.5, 49.5, 6.4, 51.51)
Coordinate Operation:
- name: Belgian Lambert 2008
- method: Lambert Conic Conformal (2SP)
Datum: European Terrestrial Reference System 1989 ensemble
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

[21]:

destination_coords = [(i.x, i.y) for i in landmarks.loc[['L001', 'L002', 'L003', 'L004'], 'geometry'].values]

Compute the transform matrix

[22]:

matrix, residual, _, _  = ggt.affine_transform_matrix(origin_coords=origin_coords, destination_coords=destination_coords)

[23]:

residual

[23]:

array([1.70078765, 5.30641372])

Save the local coordinates in wkt format in the ‘local_geometry’ column and apply the transform to the survey geodataframe

[24]:

gdf['local_geometry'] = gdf['geometry'].to_wkt()

[25]:

gdf

[25]:

	geometry	id	label	kind	class	local_geometry
0	POINT (0.000 0.000)	L001	L001	landmark	TopoPoint	POINT (0 0)
1	POINT (50.000 0.000)	L002	L002	landmark	TopoPoint	POINT (50 0)
2	POINT (4.000 120.000)	L003	L003	landmark	TopoPoint	POINT (4 120)
3	POINT (54.000 112.000)	L004	L004	landmark	TopoPoint	POINT (54 112)
4	LINESTRING (0.000 0.000, 50.000 0.000)	B01		baseline	TopoLine	LINESTRING (0 0, 50 0)
5	LINESTRING (4.000 120.000, 54.000 112.000)	B02		baseline	TopoLine	LINESTRING (4 120, 54 112)
6	POINT (2.000 0.000)			stake	TopoPoint	POINT (2 0)
7	POINT (4.000 0.000)			stake	TopoPoint	POINT (4 0)
8	POINT (6.000 0.000)			stake	TopoPoint	POINT (6 0)
9	POINT (8.000 0.000)			stake	TopoPoint	POINT (8 0)
10	POINT (10.000 0.000)			stake	TopoPoint	POINT (10 0)
11	POINT (12.000 0.000)			stake	TopoPoint	POINT (12 0)
12	POINT (14.000 0.000)			stake	TopoPoint	POINT (14 0)
13	POINT (16.000 0.000)			stake	TopoPoint	POINT (16 0)
14	POINT (18.000 0.000)			stake	TopoPoint	POINT (18 0)
15	POINT (20.000 0.000)			stake	TopoPoint	POINT (20 0)
16	POINT (22.000 0.000)			stake	TopoPoint	POINT (22 0)
17	POINT (24.000 0.000)			stake	TopoPoint	POINT (24 0)
18	POINT (26.000 0.000)			stake	TopoPoint	POINT (26 0)
19	POINT (28.000 0.000)			stake	TopoPoint	POINT (28 0)
20	POINT (30.000 0.000)			stake	TopoPoint	POINT (30 0)
21	POINT (32.000 0.000)			stake	TopoPoint	POINT (32 0)
22	POINT (34.000 0.000)			stake	TopoPoint	POINT (34 0)
23	POINT (36.000 0.000)			stake	TopoPoint	POINT (36 0)
24	POINT (5.975 119.684)			stake	TopoPoint	POINT (5.974881 119.684019)
25	POINT (7.950 119.368)			stake	TopoPoint	POINT (7.949763 119.368038)
26	POINT (9.925 119.052)			stake	TopoPoint	POINT (9.924644 119.052057)
27	POINT (11.900 118.736)			stake	TopoPoint	POINT (11.899525 118.736076)
28	POINT (13.874 118.420)			stake	TopoPoint	POINT (13.874406 118.420095)
29	POINT (15.849 118.104)			stake	TopoPoint	POINT (15.849288 118.104114)
30	POINT (17.824 117.788)			stake	TopoPoint	POINT (17.824169 117.788133)
31	POINT (19.799 117.472)			stake	TopoPoint	POINT (19.79905 117.472152)
32	POINT (21.774 117.156)			stake	TopoPoint	POINT (21.773931 117.156171)
33	POINT (23.749 116.840)			stake	TopoPoint	POINT (23.748813 116.84019)
34	POINT (25.724 116.524)			stake	TopoPoint	POINT (25.723694 116.524209)
35	POINT (27.699 116.208)			stake	TopoPoint	POINT (27.698575 116.208228)
36	POINT (29.673 115.892)			stake	TopoPoint	POINT (29.673456 115.892247)
37	POINT (31.648 115.576)			stake	TopoPoint	POINT (31.648338 115.576266)
38	POINT (33.623 115.260)			stake	TopoPoint	POINT (33.623219 115.260285)
39	POINT (35.598 114.944)			stake	TopoPoint	POINT (35.5981 114.944304)
40	POINT (37.573 114.628)			stake	TopoPoint	POINT (37.572981 114.628323)
41	POINT (39.548 114.312)			stake	TopoPoint	POINT (39.547863 114.312342)
42	LINESTRING (2.000 0.000, 5.975 119.684)	P01	P01	profile	TopoLine	LINESTRING (2 0, 5.974881 119.684019)
43	LINESTRING (6.000 0.000, 9.925 119.052)	P03	P03	profile	TopoLine	LINESTRING (6 0, 9.924644 119.052057)
44	LINESTRING (10.000 0.000, 13.874 118.420)	P05	P05	profile	TopoLine	LINESTRING (10 0, 13.874406 118.420095)
45	LINESTRING (14.000 0.000, 17.824 117.788)	P07	P07	profile	TopoLine	LINESTRING (14 0, 17.824169 117.788133)
46	LINESTRING (18.000 0.000, 21.774 117.156)	P09	P09	profile	TopoLine	LINESTRING (18 0, 21.773931 117.156171)
47	LINESTRING (22.000 0.000, 25.724 116.524)	P11	P11	profile	TopoLine	LINESTRING (22 0, 25.723694 116.524209)
48	LINESTRING (26.000 0.000, 29.673 115.892)	P13	P13	profile	TopoLine	LINESTRING (26 0, 29.673456 115.892247)
49	LINESTRING (30.000 0.000, 33.623 115.260)	P15	P15	profile	TopoLine	LINESTRING (30 0, 33.623219 115.260285)
50	LINESTRING (34.000 0.000, 37.573 114.628)	P17	P17	profile	TopoLine	LINESTRING (34 0, 37.572981 114.628323)
51	LINESTRING (7.950 119.368, 4.000 0.000)	P02	P02	profile	TopoLine	LINESTRING (7.949763 119.368038, 4 0)
52	LINESTRING (11.900 118.736, 8.000 0.000)	P04	P04	profile	TopoLine	LINESTRING (11.899525 118.736076, 8 0)
53	LINESTRING (15.849 118.104, 12.000 0.000)	P06	P06	profile	TopoLine	LINESTRING (15.849288 118.104114, 12 0)
54	LINESTRING (19.799 117.472, 16.000 0.000)	P08	P08	profile	TopoLine	LINESTRING (19.79905 117.472152, 16 0)
55	LINESTRING (23.749 116.840, 20.000 0.000)	P10	P10	profile	TopoLine	LINESTRING (23.748813 116.84019, 20 0)
56	LINESTRING (27.699 116.208, 24.000 0.000)	P12	P12	profile	TopoLine	LINESTRING (27.698575 116.208228, 24 0)
57	LINESTRING (31.648 115.576, 28.000 0.000)	P14	P14	profile	TopoLine	LINESTRING (31.648338 115.576266, 28 0)
58	LINESTRING (35.598 114.944, 32.000 0.000)	P16	P16	profile	TopoLine	LINESTRING (35.5981 114.944304, 32 0)
59	LINESTRING (39.548 114.312, 36.000 0.000)	P18	P18	profile	TopoLine	LINESTRING (39.547863 114.312342, 36 0)

[26]:

gdf['geometry'] = gpd.GeoSeries.from_wkt(gdf['local_geometry']).affine_transform(matrix)
gdf.crs = landmarks.crs

Plot the transformed survey on a basemap

NOTE:
The survey geodataframe should be converted to WGS84 (epsg:4326) for display over a basemap created with geometron.plot.plot_basemap.
It should be located entirely within the extent of the basemap

Compute the extent of the map using a buffer around the survey geodataframe

[27]:

extent = gdf.buffer(10).to_crs('epsg:4326').total_bounds[[0,2,1,3]]

[28]:

cache_dir = './tmp_files/.tiles_cache'

Plot the basemap and the survey geodataframe

[29]:

ax = plot_basemap(extent, figsize=(12,12), xyz_server='google terrain', cache_dir=cache_dir, max_zoom=19)
plot_gdf_survey(gdf.to_crs('epsg:4326'), ax=ax, extent=extent, grid='off')

[29]:

<Axes: >

../../_images/examples_geophysics_Survey_manipulate_and_display_objects_50_1.png

Export the survey geodataframe to a geopackage

[30]:

gdf.to_file('./tmp_files/survey.gpkg', driver='GPKG')