problem description¶
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, +40for Site E
profiles parallel to the main one at distace (perpendicular to the mean direction):
meters:
-20, -15, -10, -5, +5, +10, +15, +20change of CRS (done with ogr2ogr)¶
original projection:
+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defsdesired projection:
+proj=utm +zone=28 +north +ellps=GRS80 +datum=GRS80 +units=m +no_defsIn [1]:
# 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
Input and output¶
In [2]:
# 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])
In [3]:
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)
Data load¶
In [8]:
# 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'])
In [9]:
# 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)
Linear fit & perpendicular¶
In [11]:
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
Out[11]:
In [12]:
# line azimuth
# see also for a more general case
# https://github.com/marcoalopez/line_azimuth/blob/master/get_lineazi.py
Apply transform¶
In [13]:
# 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)
In [18]:
geom2shift.columns.values
Out[18]:
export data as GeoJson / shapefile / excel¶
In [19]:
# 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')
Load raster over the same area¶
In [ ]:
# TBD
Pick raster values for each point and add to geopandas¶
In [119]:
# TO BE ADDED LATER
# SEE https://github.com/mminin/snippets/blob/master/getElevationAtPoint.py
visualize materials on leaflet¶
In [20]:
import mplleaflet
In [21]:
ax = geom2shift.plot()
In [22]:
# +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)
Out[22]:
Appendix & TBD¶
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 reliefconversion 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