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()
Add the landmarks in a survey geodataframe and display the survey
[5]:
ax = plot_gdf_survey(gdf)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 1
----> 1 ax = plot_gdf_survey(gdf)
NameError: name 'gdf' is not defined
Create a baseline in local coordinates
Baseline B01 joins L001 and L002
[6]:
B01 = TopoLine(LineString([L001.geometry, L002.geometry]), id='B01', kind='baseline', show_label=False)
[7]:
B01.plot()
Baseline B02 joins L003 and L004
[8]:
B02 = TopoLine(LineString([L003.geometry, L004.geometry]), id='B02', kind='baseline', show_label=False)
Append the survey geodataframe
[9]:
gdf = pd.concat([gdf, features_to_gdf([B01, B02])], ignore_index=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[9], line 1
----> 1 gdf = pd.concat([gdf, features_to_gdf([B01, B02])], ignore_index=True)
NameError: name 'gdf' is not defined
Display the appended survey
[10]:
plot_gdf_survey(gdf)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 1
----> 1 plot_gdf_survey(gdf)
NameError: name 'gdf' is not defined
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
[11]:
number_of_profiles = 18
spacing_on_B01 = 2
spacing_on_B02 = 2
Create the stakes where the profiles start and end
[12]:
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
[13]:
gdf = pd.concat([gdf, features_to_gdf(stakes_B01), features_to_gdf(stakes_B02)], ignore_index=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 1
----> 1 gdf = pd.concat([gdf, features_to_gdf(stakes_B01), features_to_gdf(stakes_B02)], ignore_index=True)
NameError: name 'gdf' is not defined
[14]:
plot_gdf_survey(gdf)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 1
----> 1 plot_gdf_survey(gdf)
NameError: name 'gdf' is not defined
Create the profiles and append the survey geodataframe
[15]:
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'))
[16]:
gdf = pd.concat([gdf, features_to_gdf(profiles)], ignore_index=True)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 1
----> 1 gdf = pd.concat([gdf, features_to_gdf(profiles)], ignore_index=True)
NameError: name 'gdf' is not defined
[17]:
fig, ax = plt.subplots(figsize=(6,12))
plot_gdf_survey(gdf, ax=ax)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 2
1 fig, ax = plt.subplots(figsize=(6,12))
----> 2 plot_gdf_survey(gdf, ax=ax)
NameError: name 'gdf' is not defined
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
[18]:
origin_coords = [L001, L002, L003, L004]
[19]:
landmarks = gpd.read_file('./data/landmarks.shp')
landmarks.set_index('Label', drop=False, inplace=True)
landmarks.crs
[19]:
<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
[20]:
destination_coords = [(i.x, i.y) for i in landmarks.loc[['L001', 'L002', 'L003', 'L004'], 'geometry'].values]
Compute the transform matrix
[21]:
matrix, residual, _, _ = ggt.affine_transform_matrix(origin_coords=origin_coords, destination_coords=destination_coords)
[22]:
residual
[22]:
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
[23]:
gdf['local_geometry'] = gdf['geometry'].to_wkt()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 1
----> 1 gdf['local_geometry'] = gdf['geometry'].to_wkt()
NameError: name 'gdf' is not defined
[24]:
gdf
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 1
----> 1 gdf
NameError: name 'gdf' is not defined
[25]:
gdf['geometry'] = gpd.GeoSeries.from_wkt(gdf['local_geometry']).affine_transform(matrix)
gdf.crs = landmarks.crs
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[25], line 1
----> 1 gdf['geometry'] = gpd.GeoSeries.from_wkt(gdf['local_geometry']).affine_transform(matrix)
2 gdf.crs = landmarks.crs
NameError: name 'gdf' is not defined
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
[26]:
extent = gdf.buffer(10).to_crs('epsg:4326').total_bounds[[0,2,1,3]]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 1
----> 1 extent = gdf.buffer(10).to_crs('epsg:4326').total_bounds[[0,2,1,3]]
NameError: name 'gdf' is not defined
[27]:
cache_dir = './tmp_files/.tiles_cache'
Plot the basemap and the survey geodataframe
[28]:
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')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[28], line 1
----> 1 ax = plot_basemap(extent, figsize=(12,12), xyz_server='google terrain', cache_dir=cache_dir, max_zoom=19)
2 plot_gdf_survey(gdf.to_crs('epsg:4326'), ax=ax, extent=extent, grid='off')
NameError: name 'extent' is not defined
Export the survey geodataframe to a geopackage
[29]:
gdf.to_file('./tmp_files/survey.gpkg', driver='GPKG')
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[29], line 1
----> 1 gdf.to_file('./tmp_files/survey.gpkg', driver='GPKG')
NameError: name 'gdf' is not defined