%matplotlib inline
import sys,time,math
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['figure.figsize']=(8,9)
%config InlineBackend.figure_format="png"
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
from scipy.interpolate import splprep, splev
import datetime
from scipy.signal import savgol_filter
import matplotlib as mpl
import mplleaflet
from bokeh.plotting import figure, output_notebook, show
from bokeh.tile_providers import STAMEN_TERRAIN
output_notebook()
import geopandas as gp
from shapely.geometry import Point
path="/home/noel3/Dokumente/raspi-Sicherung-2018-11-25/logs/"
file_sweep="SWEEP-2018-11-25 14:58:10.717850.csv"
file_gps="GPS-2018-11-25 14:58:18.379819.csv"
dfs=pd.read_csv(path+file_sweep,header=None)
dfs.columns=["First","ANGLE","DIST","STRENGTH","timestamp"]
dfs["alpha"]=dfs.ANGLE*np.pi/180/1000.
len(dfs)
425220
dfs["time_int"]=pd.to_numeric(pd.to_datetime(dfs.timestamp))
dfs["time_int"][dfs["time_int"]<0]=np.NaN
dfs["time_fill"]=dfs[["time_int"]].interpolate(method='linear')
/home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:2: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy
dfs["time"]=pd.to_datetime(dfs.time_fill)
dfs.head()
First | ANGLE | DIST | STRENGTH | timestamp | alpha | time_int | time_fill | time | |
---|---|---|---|---|---|---|---|---|---|
0 | True | 5375 | 132 | 199 | 2018-11-25 14:58:15.964572 | 0.093811 | 1.543158e+18 | 1.543158e+18 | 2018-11-25 14:58:15.964571904 |
1 | False | 5750 | 130 | 199 | NaN | 0.100356 | NaN | 1.543158e+18 | 2018-11-25 14:58:15.967394048 |
2 | False | 6187 | 130 | 199 | NaN | 0.107984 | NaN | 1.543158e+18 | 2018-11-25 14:58:15.970216192 |
3 | False | 6562 | 130 | 199 | NaN | 0.114529 | NaN | 1.543158e+18 | 2018-11-25 14:58:15.973038336 |
4 | False | 6937 | 130 | 199 | NaN | 0.121073 | NaN | 1.543158e+18 | 2018-11-25 14:58:15.975860736 |
dfs.dtypes
First bool ANGLE int64 DIST int64 STRENGTH int64 timestamp object alpha float64 time_int float64 time_fill float64 time datetime64[ns] dtype: object
dfs=dfs[dfs.DIST>1]
dfs.index=dfs["time"]
dfs.index=dfs.index+datetime.timedelta(seconds=0.1)
dfs[80000:81010];
plt.plot(dfs[0:1000].alpha.values,"x");
#plt.plot(df[50000:50400].DIST.values/100.);
len(dfs)
389065
fig=plt.figure(figsize=(8,9));
<matplotlib.figure.Figure at 0x7f22ec9d9cc0>
def plotpolar(von,bis):
dff=dfs[von:von+bis]
ax = plt.subplot(111, projection='polar')
#ax.set_theta_direction(-1)
#ax.scatter(np.pi+dff.ANGLE.values*np.pi/180/1000.,dff.DIST.astype(float).values,c=dff.STRENGTH.values/dff.DIST.astype(float).values**2,s=1.5) #c=dff.STRENGTH.astype(float).values
ax.plot(np.pi+dff.ANGLE.values*np.pi/180/1000.,dff.DIST.astype(float).values,alpha=0.5)#,marker="o") #,c=dff.STRENGTH.values/dff.DIST.astype(float).values**2,s=1.5) #c=dff.STRENGTH.astype(float).values
#ax.set_theta_zero_location("N")
ax.set_theta_zero_location("N")
ax.set_rmax(800)
ax.grid(True)
plt.title("test"+str(dff.index[0]))
plt.show();
interact(plotpolar, von=(0,len(dfs),500),bis=(0,500,10));
coords=pd.read_csv(path+file_gps,skiprows=5)
coords.columns=["lat","lon","time"]
coords.index=pd.to_datetime(coords.time)
del coords["time"]
#coords.drop(coords.index[:1], inplace=True)
coords.head()
coords=coords[:-2]
fig=plt.figure()
plt.plot(coords.lon,coords.lat,".");
mplleaflet.display(fig=fig)
print(["coords",coords.index.min(),coords.index.max()])
print(["dfs",dfs.index.min(),dfs.index.max()])
starttime=max(coords.index.min(),dfs.index.min())
endtime=min(coords.index.max(),dfs.index.max())
cut_endtime=datetime.datetime.fromtimestamp(((starttime.value+.25*(endtime.value-starttime.value))/ 1e9-3600))
#cut_endtime=datetime.datetime.fromtimestamp(starttime.value/1e9-3600)
starttime,endtime,cut_endtime,starttime.value
['coords', Timestamp('2018-11-25 14:58:22.675404'), Timestamp('2018-11-25 15:12:52.682689')] ['dfs', Timestamp('2018-11-25 14:58:16.064571904'), Timestamp('2018-11-25 15:12:55.215343872')]
(Timestamp('2018-11-25 14:58:22.675404'), Timestamp('2018-11-25 15:12:52.682689'), datetime.datetime(2018, 11, 25, 15, 2, 0, 177225), 1543157902675404000)
dfss=pd.to_numeric(dfs[starttime:endtime].index)-starttime.value
dfss=dfss/dfss[-1]
dfs=dfs[starttime:endtime]
# spline parameters
s=0.000000030 # smoothness parameter
k=3 # spline order
nest=-1 # estimate of number of knots needed (-1 = maximal)
c=coords[starttime:endtime]
x=c.lon.values
y=c.lat.values
z=pd.to_numeric(c.index)
z=(z-z[0]).values
z=z/z[-1]
# find the knot points
tckp,u = splprep([x,y,z],s=s,k=k,nest=nest)
# evaluate spline, including interpolated points
xnew,ynew,znew = splev(dfss,tckp)
def drawmap(von,bis):
fig=plt.figure()
plt.plot(x[von:von+bis],y[von:von+bis],".");
print(von,bis)
plt.show()
interact(drawmap,von=(0,len(x)),bis=(100,500));
# derived from the Java version explained here: http://wiki.openstreetmap.org/wiki/Mercator
RADIUS = 6378137.0 # in meters on the equator
def lat2y(a):
return math.log(math.tan(math.pi / 4 + math.radians(a) / 2)) * RADIUS
def lon2x(a):
return math.radians(a) * RADIUS
lat2y = np.vectorize(lat2y)
lon2x = np.vectorize(lon2x)
coords["glat"]=coords.lat.apply(lat2y)
coords["glon"]=coords.lon.apply(lon2x)
coords.head()
lat | lon | glat | glon | |
---|---|---|---|---|
time | ||||
2018-11-25 14:58:22.675404 | 50.863256 | 12.900679 | 6.597141e+06 | 1.436097e+06 |
2018-11-25 14:58:23.674372 | 50.863259 | 12.900671 | 6.597141e+06 | 1.436096e+06 |
2018-11-25 14:58:24.674906 | 50.863258 | 12.900659 | 6.597141e+06 | 1.436095e+06 |
2018-11-25 14:58:25.675099 | 50.863253 | 12.900649 | 6.597140e+06 | 1.436094e+06 |
2018-11-25 14:58:26.675407 | 50.863247 | 12.900638 | 6.597139e+06 | 1.436092e+06 |
x_range, y_range = ((coords.glon.min(),coords.glon.max()), (coords.glat.min(),coords.glat.max()))
plot_width = int(750)
plot_height = int(plot_width//1.2)
def base_plot(tools='pan,wheel_zoom,reset',plot_width=plot_width, plot_height=plot_height, **plot_args):
p = figure(tools=tools, plot_width=plot_width, plot_height=plot_height,
x_range=x_range, y_range=y_range, outline_line_color=None,
min_border=0, min_border_left=0, min_border_right=0,
min_border_top=0, min_border_bottom=0, **plot_args)
p.axis.visible = False
p.xgrid.grid_line_color = None
p.ygrid.grid_line_color = None
return p
options = dict(line_color=None, fill_color='blue', size=5)
orient=-np.arctan2(savgol_filter(lon2x(xnew),225,1,deriv=1),savgol_filter(lat2y(ynew),225,1,deriv=1))+np.pi/2
plt.plot(orient)
plt.plot((xnew-xnew[0])/xnew.std())
plt.plot((ynew-ynew[0])/ynew.std());
fig=plt.figure()
fit,=plt.plot(xnew[::200],ynew[::200],'xr',label='fit')
plt.quiver(xnew[0::2000],ynew[0::2000],np.cos(orient[0::2000]),np.sin(orient[0::2000]), units='width',scale=20.2,pivot="middle");
plt.axes().set_aspect('equal', 'datalim')
mplleaflet.display(fig=fig)
dfs.head()
First | ANGLE | DIST | STRENGTH | timestamp | alpha | time_int | time_fill | time | |
---|---|---|---|---|---|---|---|---|---|
time | |||||||||
2018-11-25 14:58:22.681509120 | False | 181000 | 3852 | 20 | NaN | 3.159046 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.581509120 |
2018-11-25 14:58:22.704183040 | False | 205812 | 5352 | 25 | NaN | 3.592097 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.604183040 |
2018-11-25 14:58:22.726856704 | False | 230375 | 3972 | 25 | NaN | 4.020802 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.626856704 |
2018-11-25 14:58:22.749530368 | False | 255062 | 4293 | 25 | NaN | 4.451672 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.649530368 |
2018-11-25 14:58:22.772204032 | False | 279750 | 7587 | 25 | NaN | 4.882559 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.672204032 |
k=np.pi/180
k2=360/(100.*RADIUS*2*np.pi)
Scanseorient=-np.pi/2
#y-Drehung
x1=0
y1=dfs.DIST.values*np.cos(Scanseorient+dfs.alpha.values)
z1=dfs.DIST.values*np.sin(Scanseorient+dfs.alpha.values)
#z-Drehung
x2=x1*np.cos(orient)-y1*np.sin(orient)
y2=x1*np.sin(orient)+y1*np.cos(orient)
z2=z1
#loc to global
##xtar=lon2x(xnew+k2*dfs.DIST.values*np.cos(orient+dfs.alpha.values)/np.cos(ynew[0]))
##ytar=lat2y(ynew+k2*dfs.DIST.values*np.sin(orient+dfs.alpha.values))
xtar=lon2x(xnew+k2*x2/np.cos(ynew[0]))
ytar=lat2y(ynew+k2*y2)
ztar=np.clip(z2/100,-3,1000)
von,bis=00000,60000
p = base_plot()
p.add_tile(STAMEN_TERRAIN)
p.circle(x=xtar[von:bis:5], y=ytar[von:bis:5],alpha=.09, **options)
show(p)
pd.to_numeric(dfs[von:bis].index).values
array([1543158324457358336, 1543158324459362560, 1543158324461366528, ..., 1543158428279106816, 1543158428281395456, 1543158428283683840])
plt.plot(dfs.index,np.linspace(0,1,len(dfs)),label="Sweep")
plt.plot(coords.index,np.linspace(0,1,len(coords)),label="GPS")
plt.legend();
path_po="/home/noel3/dev/workspaces/PotreeConverter/master/test/resources_moved/"
file_po="ripple_ascii.ply"
df=pd.read_csv(path+file,skiprows=15,sep=" ",header=None)
df.columns=["x","y","z","r","g","b","nx","ny","nz","i"]
df.x=0.00001*df.x+12.88370
df.y=0.00001*df.y+50.86434
df.z=1*df.z
df.head()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-149-54fc39f5b10c> in <module>() ----> 1 df=pd.read_csv(path+file,skiprows=15,sep=" ",header=None) 2 df.columns=["x","y","z","r","g","b","nx","ny","nz","i"] 3 4 df.x=0.00001*df.x+12.88370 5 df.y=0.00001*df.y+50.86434 NameError: name 'file' is not defined
dfs["x"]=xtar
dfs["y"]=ytar
dfs["z"]=ztar
dfs.head()
First | ANGLE | DIST | STRENGTH | timestamp | alpha | time_int | time_fill | time | y | x | z | Coordinates | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
time | |||||||||||||
2018-11-25 14:58:22.681509120 | False | 181000 | 3852 | 20 | NaN | 3.159046 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.581509120 | 6.597142e+06 | 1.436097e+06 | 38.514133 | POINT (1436097.054111757 6597142.333480542) |
2018-11-25 14:58:22.704183040 | False | 205812 | 5352 | 25 | NaN | 3.592097 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.604183040 | 6.597178e+06 | 1.436092e+06 | 48.180181 | POINT (1436091.810708485 6597177.528666587) |
2018-11-25 14:58:22.726856704 | False | 230375 | 3972 | 25 | NaN | 4.020802 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.626856704 | 6.597189e+06 | 1.436090e+06 | 25.331832 | POINT (1436090.099939776 6597188.862096393) |
2018-11-25 14:58:22.749530368 | False | 255062 | 4293 | 25 | NaN | 4.451672 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.649530368 | 6.597206e+06 | 1.436088e+06 | 11.066223 | POINT (1436087.561378659 6597205.787769186) |
2018-11-25 14:58:22.772204032 | False | 279750 | 7587 | 25 | NaN | 4.882559 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.672204032 | 6.597258e+06 | 1.436080e+06 | -3.000000 | POINT (1436079.862828478 6597257.569312509) |
dfs.columns
Index(['First', 'ANGLE', 'DIST', 'STRENGTH', 'timestamp', 'alpha', 'time_int', 'time_fill', 'time', 'y', 'x', 'z', 'Coordinates'], dtype='object')
dfs['Coordinates'] = list(zip(dfs.x, dfs.y))
dfs['Coordinates'] = dfs['Coordinates'].apply(Point)
gdf = gp.GeoDataFrame(dfs, geometry='Coordinates')
#gdf.plot(alpha=.1);
gdf.crs = {'init' :'epsg:3857'}
gdf.crs
{'init': 'epsg:3857'}
gdf2=gdf.to_crs({'init': 'epsg:5243'}) #https://epsg.io/5243
gdf2['x'] = gdf2['Coordinates'].apply(lambda x: x.coords[0][0])
gdf2['y'] = gdf2['Coordinates'].apply(lambda x: x.coords[0][1])
gdf2=gdf2.drop(labels="Coordinates",axis=1)
gdf2.head()
First | ANGLE | DIST | STRENGTH | timestamp | alpha | time_int | time_fill | time | y | x | z | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
time | ||||||||||||
2018-11-25 14:58:22.681509120 | False | 181000 | 3852 | 20 | NaN | 3.159046 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.581509120 | -12440.808321 | 168829.184460 | 38.514133 |
2018-11-25 14:58:22.704183040 | False | 205812 | 5352 | 25 | NaN | 3.592097 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.604183040 | -12418.749533 | 168825.149163 | 48.180181 |
2018-11-25 14:58:22.726856704 | False | 230375 | 3972 | 25 | NaN | 4.020802 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.626856704 | -12411.646720 | 168823.835651 | 25.331832 |
2018-11-25 14:58:22.749530368 | False | 255062 | 4293 | 25 | NaN | 4.451672 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.649530368 | -12401.038852 | 168821.884345 | 11.066223 |
2018-11-25 14:58:22.772204032 | False | 279750 | 7587 | 25 | NaN | 4.882559 | NaN | 1.543158e+18 | 2018-11-25 14:58:22.672204032 | -12368.584416 | 168815.957478 | -3.000000 |
#x y z r g b nx ny nz i
gdf3=gdf2[["x","y","z"]]
#gdf3["z"]=0.0015*gdf2["DIST"]*np.abs(np.sin(gdf2["ANGLE"]/(2**16)))
gdf3["r"]=gdf2["STRENGTH"]%50*5
gdf3["g"]=0
gdf3["b"]=255-gdf2["STRENGTH"]%50*5
gdf3["nx"]=0.3;gdf3["ny"]=0.3;gdf3["nz"]=0.3
gdf3["i"]=gdf2["STRENGTH"]
gdf3.head()
out=gdf3.to_csv(header=None,index=None,sep=" ",float_format='%.6f')
/home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:5: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy """ /home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:6: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy /home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:7: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy import sys /home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:9: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy if __name__ == '__main__': /home/noel3/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:10: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/indexing.html#indexing-view-versus-copy # Remove the CWD from sys.path while we load stuff.
out_start='''\
ply
format ascii 1.0
comment Author: CloudCompare (TELECOM PARISTECH/EDF R&D) moved by T Alder
obj_info Generated by CloudCompare!
element vertex '''+str(len(gdf3))+'''
property float x
property float y
property float z
property uchar red
property uchar green
property uchar blue
property float nx
property float ny
property float nz
end_header
'''
f = open(path_po+"ripple_ascii_tim-vertical02.ply", "w")
f.write(out_start+out)
29969070
!~/dev/workspaces/PotreeConverter/master/build/PotreeConverter/PotreeConverter \
"/home/noel3/dev/workspaces/PotreeConverter/master/test/resources_moved/ripple_ascii_tim-vertical02.ply" \
-o "/home/noel3/dev/workspaces/PotreeConverter/master/test/converted52" \
-p test51 -projection "+proj=lcc +lat_1=48.66666666666666 +lat_2=53.66666666666666 +lat_0=51 +lon_0=10.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"\
--overwrite -l 5 --output-format LAZ
== params == source[0]: /home/noel3/dev/workspaces/PotreeConverter/master/test/resources_moved/ripple_ascii_tim-vertical02.ply outdir: /home/noel3/dev/workspaces/PotreeConverter/master/test/converted52 spacing: 0 diagonal-fraction: 200 levels: 5 format: scale: 0 pageName: test51 projection: +proj=lcc +lat_1=48.66666666666666 +lat_2=53.66666666666666 +lat_0=51 +lon_0=10.5 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs--overwrite AABB: min: [167.574, -13.228,2, -3] max: [168.893, -12.368,6, 100,523] size: [1.319,7, 859,575, 103,523] cubic AABB: min: [167.574, -13.228,2, -3] max: [168.893, -11.908,5, 1.316,7] size: [1.319,7, 1.319,7, 1.319,7] spacing calculated from diagonal: 11,429 READING: /home/noel3/dev/workspaces/PotreeConverter/master/test/resources_moved/ripple_ascii_tim-vertical02.ply closing writer conversion finished 385.106 points were processed and 385.106 points ( 100% ) were written to the output. duration: 4,575s
#["#%02x%02x%02x" % (int(r), int(g), int(b)) for r, g, b, _ in 255*mpl.cm.jet(mpl.colors.Normalize()(dfs[von:bis].STRENGTH.values/dfs[von:bis].DIST.values))]
(0.015*gdf2["DIST"][:1000]*np.cos(gdf2["ANGLE"][:1000]/(2**16))).plot();
#((gdf2["ANGLE"][:1000]/(2**16))).plot(secondary_y=True);