The approach might be useful for those working on data collected over sub-linear profiles, e.g. for planning or implementation purposes...
for Site D:
profiles parallel to the main one at distace (perpendicular to the mean direction):
meters:
-40, -30, -20, -10, +10, +20, +30, +40
for Site E
profiles parallel to the main one at distace (perpendicular to the mean direction):
meters:
-20, -15, -10, -5, +5, +10, +15, +20
original projection:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
desired projection:
+proj=utm +zone=28 +north +ellps=GRS80 +datum=GRS80 +units=m +no_defs
# files as WGS84 geographic, as exported from RTK processing
input_raw_D = './0_location-profiles-RAW/AGPA-G_D2.shp'
input_raw_E = './0_location-profiles-RAW/AGPA-E.shp'
output_UTM_D = './1_location-profiles-UTM/ERT-siteD-UTM28-GRS80.shp'
# get_input projection
import fiona
from fiona.crs import to_string
c = fiona.open(output_UTM_D, 'r')
# print(c.crs)
print(to_string(c.crs))
# wgs = from_epsg
# wgs_proj4_string = to_string(wgs)
# print(wgs_proj4_string)
UTM_D = './0_location-profiles-RAW/AGPA-G_D2.shp'
UTM_E = './0_location-profiles-RAW/AGPA-E.shp'
# one could reproject via
# https://pcjericks.github.io/py-gdalogr-cookbook/projection.html#reproject-a-geometry
+ellps=GRS80 +no_defs +proj=utm +units=m +zone=28
# arrays for shifted points
delta_D = [-40, -30, -20, -10, +10, +20, +30, +40]
delta_E = [-20, -15, -10, -5, +5, +10, +15, +20]
delta_TEST = [10]
print(delta_D)
print(delta_E)
print(delta_D[0])
[-40, -30, -20, -10, 10, 20, 30, 40] [-20, -15, -10, -5, 5, 10, 15, 20] -40
delta_select = 'D'
delta = delta_D
# input
input_file = './1_location-profiles-UTM/ERT-site' + delta_select + '-UTM28-GRS80.shp'
# output CSV/json
outcsv = "./2_output/agpa_" + delta_select + "-shifted.csv"
outjson = './2_output/agpa_' + delta_select + '-shifted.geojson'
# delta
# outjson = '2_output/test.geojson'
# outjson = '2_output/agpa_D-shifted.geojson'
# outjson = '2_output/agpa_E-shifted.geojson'
print(delta_select)
print(input_file)
print(outcsv)
print(outjson)
print(delta)
D ./1_location-profiles-UTM/ERT-siteD-UTM28-GRS80.shp ./2_output/agpa_D-shifted.csv ./2_output/agpa_D-shifted.geojson [-40, -30, -20, -10, 10, 20, 30, 40]
# load data
import fiona
import shapely
c = fiona.open(input_file, 'r')
c.crs
# number of points, records
len(c) # should be 48 for AGPA-G
# PRINT SCHEMA
print('schema is:')
print(c.schema, '\n')
# PRINT GEOMETRIES
print('geometries are: ')
for feat in c:
print(feat['geometry'])
schema is: {'properties': OrderedDict([('antenna he', 'float:24.15'), ('lateral rm', 'float:24.15'), ('name', 'str:80'), ('projected', 'str:80'), ('rms', 'str:80'), ('sample cou', 'float:24.15'), ('solution s', 'str:80')]), 'geometry': 'Point'} geometries are: {'type': 'Point', 'coordinates': (648601.4745557815, 3229147.9738091887)} {'type': 'Point', 'coordinates': (648600.5361209628, 3229143.138978075)} {'type': 'Point', 'coordinates': (648599.7661030253, 3229138.156342755)} {'type': 'Point', 'coordinates': (648598.7990097381, 3229133.272981949)} {'type': 'Point', 'coordinates': (648597.9483836357, 3229128.388791427)} {'type': 'Point', 'coordinates': (648596.9545805457, 3229123.46105742)} {'type': 'Point', 'coordinates': (648596.0947405221, 3229118.56331744)} {'type': 'Point', 'coordinates': (648595.1762005028, 3229113.609735632)} {'type': 'Point', 'coordinates': (648594.2721767295, 3229108.636121245)} {'type': 'Point', 'coordinates': (648593.2952063341, 3229103.763752384)} {'type': 'Point', 'coordinates': (648592.6392054026, 3229098.807557541)} {'type': 'Point', 'coordinates': (648591.786515621, 3229094.036386801)} {'type': 'Point', 'coordinates': (648590.9748514218, 3229088.9411467803)} {'type': 'Point', 'coordinates': (648589.9975829801, 3229084.1011719736)} {'type': 'Point', 'coordinates': (648589.144341087, 3229079.2450881554)} {'type': 'Point', 'coordinates': (648588.2456953059, 3229074.411880112)} {'type': 'Point', 'coordinates': (648587.5794628204, 3229069.5876915418)} {'type': 'Point', 'coordinates': (648586.5790606767, 3229064.5499946526)} {'type': 'Point', 'coordinates': (648585.8623549817, 3229059.698842795)} {'type': 'Point', 'coordinates': (648585.0456319387, 3229054.725776453)} {'type': 'Point', 'coordinates': (648584.3042273485, 3229050.012683716)} {'type': 'Point', 'coordinates': (648583.4001623433, 3229045.0961384033)} {'type': 'Point', 'coordinates': (648582.5647751306, 3229039.665844577)} {'type': 'Point', 'coordinates': (648581.665703492, 3229035.211213602)} {'type': 'Point', 'coordinates': (648580.9057197726, 3229030.0542368554)} {'type': 'Point', 'coordinates': (648579.9821064426, 3229025.2335872888)} {'type': 'Point', 'coordinates': (648579.2077628493, 3229020.400097237)} {'type': 'Point', 'coordinates': (648578.3236284078, 3229015.495011512)} {'type': 'Point', 'coordinates': (648577.4638452544, 3229010.754280825)} {'type': 'Point', 'coordinates': (648576.5390187646, 3229005.727932812)} {'type': 'Point', 'coordinates': (648575.5946259949, 3229000.8785012243)} {'type': 'Point', 'coordinates': (648574.9364234685, 3228995.664924923)} {'type': 'Point', 'coordinates': (648574.0533106204, 3228991.0935741314)} {'type': 'Point', 'coordinates': (648573.1830404663, 3228986.0321496315)} {'type': 'Point', 'coordinates': (648572.2467538107, 3228980.9361673873)} {'type': 'Point', 'coordinates': (648571.4503726437, 3228976.1057025185)} {'type': 'Point', 'coordinates': (648570.5414256933, 3228970.8126438186)} {'type': 'Point', 'coordinates': (648569.672187713, 3228966.050441338)} {'type': 'Point', 'coordinates': (648568.3397000737, 3228961.1769017302)} {'type': 'Point', 'coordinates': (648567.7542779779, 3228956.3556339275)} {'type': 'Point', 'coordinates': (648566.9048521187, 3228951.458135281)} {'type': 'Point', 'coordinates': (648566.109367921, 3228946.6162900357)} {'type': 'Point', 'coordinates': (648565.0290407431, 3228941.8700283184)} {'type': 'Point', 'coordinates': (648563.3972566999, 3228936.4794101315)} {'type': 'Point', 'coordinates': (648562.8736383251, 3228931.9387767552)} {'type': 'Point', 'coordinates': (648561.503112745, 3228927.599268888)} {'type': 'Point', 'coordinates': (648560.9811715708, 3228922.554788729)} {'type': 'Point', 'coordinates': (648559.6440178016, 3228917.816313375)}
# IMPORT GEOMETRIES IN GEOPANDAS
import pandas as pd
import geopandas as gpd
gdf = gpd.GeoDataFrame.from_file(input_file)
# gdf.head()
# create empty dataframe
columns = ['x_UTM','y_UTM']
dictionary_append = []
geom2fit = pd.DataFrame(columns=columns)
# print all rows x,y
# fill a dataframe for each row
for geometry in gdf.geometry:
print(geometry.x)
print(geometry.y)
dictionary_append.append({'x_UTM':geometry.x, 'y_UTM':geometry.y})
geom2fit = geom2fit.append(dictionary_append)
print(geom2fit.x_UTM.values, geom2fit.y_UTM.values)
648601.4745557815 3229147.9738091887 648600.5361209628 3229143.138978075 648599.7661030253 3229138.156342755 648598.7990097381 3229133.272981949 648597.9483836357 3229128.388791427 648596.9545805457 3229123.46105742 648596.0947405221 3229118.56331744 648595.1762005028 3229113.609735632 648594.2721767295 3229108.636121245 648593.2952063341 3229103.763752384 648592.6392054026 3229098.807557541 648591.786515621 3229094.036386801 648590.9748514218 3229088.9411467803 648589.9975829801 3229084.1011719736 648589.144341087 3229079.2450881554 648588.2456953059 3229074.411880112 648587.5794628204 3229069.5876915418 648586.5790606767 3229064.5499946526 648585.8623549817 3229059.698842795 648585.0456319387 3229054.725776453 648584.3042273485 3229050.012683716 648583.4001623433 3229045.0961384033 648582.5647751306 3229039.665844577 648581.665703492 3229035.211213602 648580.9057197726 3229030.0542368554 648579.9821064426 3229025.2335872888 648579.2077628493 3229020.400097237 648578.3236284078 3229015.495011512 648577.4638452544 3229010.754280825 648576.5390187646 3229005.727932812 648575.5946259949 3229000.8785012243 648574.9364234685 3228995.664924923 648574.0533106204 3228991.0935741314 648573.1830404663 3228986.0321496315 648572.2467538107 3228980.9361673873 648571.4503726437 3228976.1057025185 648570.5414256933 3228970.8126438186 648569.672187713 3228966.050441338 648568.3397000737 3228961.1769017302 648567.7542779779 3228956.3556339275 648566.9048521187 3228951.458135281 648566.109367921 3228946.6162900357 648565.0290407431 3228941.8700283184 648563.3972566999 3228936.4794101315 648562.8736383251 3228931.9387767552 648561.503112745 3228927.599268888 648560.9811715708 3228922.554788729 648559.6440178016 3228917.816313375 [648601.47455578 648600.53612096 648599.76610303 648598.79900974 648597.94838364 648596.95458055 648596.09474052 648595.1762005 648594.27217673 648593.29520633 648592.6392054 648591.78651562 648590.97485142 648589.99758298 648589.14434109 648588.24569531 648587.57946282 648586.57906068 648585.86235498 648585.04563194 648584.30422735 648583.40016234 648582.56477513 648581.66570349 648580.90571977 648579.98210644 648579.20776285 648578.32362841 648577.46384525 648576.53901876 648575.59462599 648574.93642347 648574.05331062 648573.18304047 648572.24675381 648571.45037264 648570.54142569 648569.67218771 648568.33970007 648567.75427798 648566.90485212 648566.10936792 648565.02904074 648563.3972567 648562.87363833 648561.50311275 648560.98117157 648559.6440178 ] [3229147.97380919 3229143.13897807 3229138.15634276 3229133.27298195 3229128.38879143 3229123.46105742 3229118.56331744 3229113.60973563 3229108.63612125 3229103.76375238 3229098.80755754 3229094.0363868 3229088.94114678 3229084.10117197 3229079.24508816 3229074.41188011 3229069.58769154 3229064.54999465 3229059.6988428 3229054.72577645 3229050.01268372 3229045.0961384 3229039.66584458 3229035.2112136 3229030.05423686 3229025.23358729 3229020.40009724 3229015.49501151 3229010.75428083 3229005.72793281 3229000.87850122 3228995.66492492 3228991.09357413 3228986.03214963 3228980.93616739 3228976.10570252 3228970.81264382 3228966.05044134 3228961.17690173 3228956.35563393 3228951.45813528 3228946.61629004 3228941.87002832 3228936.47941013 3228931.93877676 3228927.59926889 3228922.55478873 3228917.81631337]
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from math import atan, degrees
# x of geometry in UTM
x_points = geom2fit.x_UTM.values*1.0
# y of geometry in UTM
y_points = geom2fit.y_UTM.values*1.0
print(x_points)
# center point
center_point = [np.mean(x_points), np.mean(y_points)]
print(center_point)
# center_point_x = np.mean(x_points)
center_point_x = (max(x_points) + min(x_points))/2
# center_point_y = np.mean(y_points)
center_point_y = (max(y_points) + min(y_points))/2
print('center point x,y are :', center_point_x, center_point_y)
# fit of points
fit = np.polyfit(x_points,y_points,1)
m, b = np.polyfit(x_points, y_points, 1)
fit_fn = np.poly1d(fit)
# print(fit)
print('slope of fit is: ', fit[0])
print('intercept of fit is: ', fit[1])
### perpendicular line passing through the center
print('slope of perpendicular to fit is: ', -1/fit[0])
m_p = (-1/fit[0])
# m_angle = degrees(atan(fit[0]))
# m_p_angle = degrees(atan(m_p))
m_angle = atan(fit[0])
m_p_angle = atan(m_p)
print('slope is: ', m_angle)
print('perpendicular slope is: ', m_p_angle) #in radians
print('check: ', m_angle - m_p_angle) # in radians
# intercept of perpendicular line
b_c = center_point_y -m_p*center_point_x
y_p_points = m_p*x_points + b_c
# set figure sixe
fig = plt.gcf()
fig.set_size_inches(4.5, 10.5)
# perpendicular_fit = ((-1/m)*x_points, y_points)
# points and fit
ax1 = plt.plot(x_points,y_points, 'yo', x_points, fit_fn(x_points), '--k')
# plot center point
ax2 = plt.plot(center_point_x, center_point_y, 'o')
# plot perpendicular line
ax3 = plt.plot(x_points, y_p_points)
# plot perpendicular fit
# plt.plot(x_points,(-1/m)*x_points, 'ro')
plt.grid('on')
plt.axes().set_aspect(1) # it's turned off, but lines are perpendicular, uncomment to check ;)
plt.show
[648601.47455578 648600.53612096 648599.76610303 648598.79900974 648597.94838364 648596.95458055 648596.09474052 648595.1762005 648594.27217673 648593.29520633 648592.6392054 648591.78651562 648590.97485142 648589.99758298 648589.14434109 648588.24569531 648587.57946282 648586.57906068 648585.86235498 648585.04563194 648584.30422735 648583.40016234 648582.56477513 648581.66570349 648580.90571977 648579.98210644 648579.20776285 648578.32362841 648577.46384525 648576.53901876 648575.59462599 648574.93642347 648574.05331062 648573.18304047 648572.24675381 648571.45037264 648570.54142569 648569.67218771 648568.33970007 648567.75427798 648566.90485212 648566.10936792 648565.02904074 648563.3972567 648562.87363833 648561.50311275 648560.98117157 648559.6440178 ] [648581.0571522126, 3229032.6283563185] center point x,y are : 648580.5592867916 3229032.8950612815 slope of fit is: 5.643287879511387 intercept of fit is: -431096.99035144487 slope of perpendicular to fit is: -0.17720166352501993 slope is: 1.3954152207670156 perpendicular slope is: -0.17538110602788104 check: 1.5707963267948966
/Users/anrossi/miniconda/envs/python3test/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance. warnings.warn(message, mplDeprecation, stacklevel=1)
<function matplotlib.pyplot.show>
# line azimuth
# see also for a more general case
# https://github.com/marcoalopez/line_azimuth/blob/master/get_lineazi.py
# translate the center point according to the delta_D
from math import sqrt, sin, cos, tan
# from geopandas import GeoDataFrame
import pandas as pd
import geopandas as gpd
# from shapely.geometry import Point
import shapely
#####
# FROM ABOVE:
### perpendicular line passing through the center
print('slope of perpendicular to fit is: ', -1/fit[0])
m_p = (-1/fit[0])
# intercept of perpendicular line
b_c = center_point_y -m_p*center_point_x
y_p_points = m_p*x_points + b_c
#####
# test, put points at given distance along the perpendicular line
# delta_D = delta_TEST
# create a dataframe
# columns = ['x_UTM','y_UTM']
# dictionary_append = []
# geom2shift = pd.DataFrame(columns=columns)
'''
# create a dataframe
columns = ['x_UTM','y_UTM','sPoints','delta_D']
dictionary_append = []
#geom2shift = pd.DataFrame(columns=columns)
#geom2shift.set_geometry('sPoints')
geom2shift = gpd.GeoDataFrame(columns=columns,geometry='sPoints')
print('delta_D is: ',delta_D)
'''
# fix by Mikhail Minin start
# to define the data frame
columns = ['x_UTM','y_UTM','sPoints','delta']
dictionary_append = []
#geom2shift = pd.DataFrame(columns=columns)
#geom2shift.set_geometry('sPoints')
geom2shift = gpd.GeoDataFrame(columns=columns,geometry='sPoints')
# fix by Mikhail Minin stop
# delta = delta_E
print('delta is: ',delta)
j=0
for element in delta:
print(j)
print(element)
for i in range(len(x_points)):
y_shift = 0
x_shift = 0
x_shift = (element * cos(m_p_angle))
y_shift = (element * sin(m_p_angle))
# print(x_points[0])
# print(y_points[0])
# print('x_shift is: ',x_shift)
# print('y_shift is: ',y_shift)
x_shifted = x_points + x_shift
y_shifted = y_points + y_shift
# print('x_shifted is: ',x_shifted)
# print('y_shifted is: ',y_shifted)
geom2shift.loc[j+i,'sPoints']=shapely.geometry.point.Point(x_shifted[i],y_shifted[i])
geom2shift.loc[j+i,'x_UTM']=x_shifted[i]
geom2shift.loc[j+i,'y_UTM']=y_shifted[i]
geom2shift.loc[j+i,'delta']=element
j+=len(x_points)
# calculate x, y for delta_d elements
# plot each point shifted
# points and fit
fig.set_size_inches(4.5, 10.5)
ax1 = plt.plot(x_points, fit_fn(x_points), '--k')
fig = plt.gcf()
ax2 = plt.plot(x_points[24],y_points[24])
plt.plot(x_shifted,y_shifted, '.', color='blue')
# plot perpendicular line
ax3 = plt.plot(x_points, y_points, '.', color='red')
plt.grid('on')
plt.axes().set_aspect(1) # it's turned off, but lines are perpendicular, uncomment to check ;)
# plot perpendicular line
ax3 = plt.plot(x_points, y_p_points)
# for i in range(len(x_points)):
# geom2shift.loc[j+i,'sPoints']=shapely.geometry.point.Point(x_points[i],y_p_points[i])
# geom2shift.loc[j+i,'x_UTM']=x_points[i]
# geom2shift.loc[j+i,'y_UTM']=y_p_points[i]
# geom2shift.loc[j+i,'delta_D']=element
# j+=len(x_points)
plt.show
# populate the pandas dataframe with new results
# for geometry in geom2shift.geometry:
# print(geometry)
# add geometry
# add attributes based on element
# save
# # create empty dataframe
# columns = ['x_UTM','y_UTM']
# dictionary_append = []
# geom2fit = pd.DataFrame(columns=columns)
# # print all rows x,y
# # fill a dataframe for each row
# for geometry in gdf.geometry:
# print(geometry.x)
# print(geometry.y)
# dictionary_append.append({'x_UTM':geometry.x, 'y_UTM':geometry.y})
# geom2fit = geom2fit.append(dictionary_append)
slope of perpendicular to fit is: -0.17720166352501993 delta is: [-40, -30, -20, -10, 10, 20, 30, 40] 0 -40 48 -30 96 -20
/Users/anrossi/miniconda/envs/python3test/lib/python3.6/site-packages/matplotlib/cbook/deprecation.py:106: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance. In a future version, a new instance will always be created and returned. Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance. warnings.warn(message, mplDeprecation, stacklevel=1)
144 -10 192 10 240 20 288 30 336 40
geom2shift.columns.values
array(['x_UTM', 'y_UTM', 'sPoints', 'delta'], dtype=object)
# export geopandas as json/shapefile
# outjson = '2_output/test.geojson'
# outjson = '2_output/agpa_D-shifted.geojson'
# outjson = '2_output/agpa_E-shifted.geojson'
with open(outjson, 'w') as f:
f.write(geom2shift.to_json())
# export to csv
with open(outcsv, 'w') as f:
f.write(geom2shift.to_csv())
# df.to_file('/dev/sdc/Boundary.shp', driver='ESRI Shapefile', crs_wkt=prj)
# and to csv as
# geom2shift.to_csv()
# to convert to shapefile you'd need an appropriate driver installed, but it should look something like
# geom2shift.to_file(driver ='ESRI Shapefile', filename='abc.shp')
# TBD
# TO BE ADDED LATER
# SEE https://github.com/mminin/snippets/blob/master/getElevationAtPoint.py
import mplleaflet
ax = geom2shift.plot()
# +ellps=GRS80 +no_defs +proj=utm +units=m +zone=28
crs_leaflet = {'ellps': 'GRS80',
'proj': 'utm',
'units': 'm',
'zone': '28'}
mplleaflet.display(fig=ax.figure, crs=crs_leaflet)
raster conversion to UTM GRS80
# dtm
gdalwarp -t_srs "+proj=utm +zone=28 +north +ellps=GRS80 \
+datum=GRS80 +units=m +no_defs" \
dtm-lanzarote-wgs84-utm28n.tif dtm-lanzarote-grs80-utm28n.tif
# same for shade relief
conversion of geojson to shapefile (table to contain all info on point locations)
ogr2ogr -nlt POINT -skipfailures D.shp agpa_D-shifted.geojson OGRGeoJSON
ogr2ogr -nlt POINT -skipfailures E.shp agpa_E-shifted.geojson OGRGeoJSON