# -*- coding: utf-8 -*- """ Created on Wed Oct 6 14:39:58 2021 @author: karaouli """ import matplotlib.pyplot as plt import geopandas as gpd import numpy as np from matplotlib import path import fiona from shapely import geometry import pyproj import sqlite3 import pandas as pd from tqdm import tqdm import pickle import sys sys.path.insert(0, r'D:\karaouli\Desktop\Projects\dutch_dabases\module\geo_vtk\src') from vtkclass import VtkClass import cpt_module class Bro_Reader: def __init__(self,datapath=r'.\data\\cpts\\',polygon=None,Xmin=None,Xmax=None,Ymin=None,Ymax=None): """ This package reads the gpkg file with cpts (from BRO), makes a selection based on the polygon provided, and plots all CPTs in VTK format (to be viewed with Paraview) All files will be stored to the results\data for the selected file (pickle fomrat) and results\vtk\cpt.vtk for the graphics output. Parameters ---------- :datapath Provide path to the gpkg file :polygon Fiona shapefile. :Xmin :Xmax :Ymin :Ymax Returns ------- None. """ self.datapath=datapath self.cpts=[] self.wgs84_p = pyproj.Proj(init='epsg:4326') self.utm_xy_31 = pyproj.Proj(init='epsg:32631') self.syn_21 = pyproj.Proj(init='epsg:3414') self.rd_p = pyproj.Proj(init='epsg:28992') self.results=[] if Xmin==None: # make poygon from coordinates # self.p = path.Path([(-np.Inf,-np.Inf), (-np.Inf, np.Inf), (np.Inf, np.Inf), (-np.Inf, -np.Inf)]) # square with legs length 1 and bottom left corner at the origin p1 = geometry.Point(-np.Inf,-np.Inf) p2 = geometry.Point(np.Inf,-np.Inf) p3 = geometry.Point(np.Inf,np.Inf) p4 = geometry.Point(-np.Inf,np.Inf) pointList = [p1, p2, p3, p4, p1] self.p= geometry.Polygon([[p.x, p.y] for p in pointList]) else: self.Xmin=Xmin self.Xmax=Xmax self.Ymin=Ymin self.Ymax=Ymax self.p=self.make_bbox() def make_bbox(self): # convert to lat,lon data=np.c_[[self.Xmin,self.Xmax],[self.Ymin,self.Ymax]] ulx,uly=pyproj.transform(self.rd_p,self.wgs84_p, data[:,0], data[:,1]) Xmin=np.min(ulx) Xmax=np.max(ulx) Ymin=np.min(uly) Ymax=np.max(uly) p1 = geometry.Point(Xmin,Ymin) p2 = geometry.Point(Xmax,Ymin) p3 = geometry.Point(Xmax,Ymax) p4 = geometry.Point(Xmin,Ymax) pointList = [p1, p2, p3, p4, p1] p= geometry.Polygon([[p.x, p.y] for p in pointList]) return p def read_cpt(self): # attribures=fiona.listlayers(self.datapath+"brocptvolledigeset.gpkg") cpts=gpd.read_file(self.datapath+"brocptvolledigeset.gpkg",mask=self.p) # cpts=gpd.read_file(self.datapath+"brocptvolledigeset.gpkg",mask=self.p,layer=['cpt_cone_penetration_test_result']['depth']) # flags = self.p.contains_points(np.c_[cpts['x_or_lon'].to_numpy(dtype='float'),cpts['y_or_lat'].to_numpy(dtype='float')]) # plt.plot(*self.p.exterior.xy) # self.cpts=cpts.iloc[flags] conn = sqlite3.connect(self.datapath+"brocptvolledigeset.gpkg") cursor = conn.cursor() sql="SELECT cpt_cone_penetration_test_result.penetration_length,\ cpt_cone_penetration_test_result.cone_resistance,\ cpt_cone_penetration_test_result.depth,\ cpt_cone_penetration_test_result.local_friction,\ cpt_cone_penetration_test_result.friction_ratio,\ cpt_cone_penetration_test_result.inclination_resultant,\ cpt_cone_penetration_test_result.pore_pressure_u1,\ cpt_cone_penetration_test_result.pore_pressure_u2,\ cpt_cone_penetration_test_result.pore_pressure_u3,\ cpt_geotechnical_survey.x_or_lon,\ cpt_geotechnical_survey.y_or_lat,\ cpt_geotechnical_survey.offset,\ cpt_geotechnical_survey.bro_id,\ cpt_cone_penetrometer.cone_surface_quotient\ FROM cpt_geotechnical_survey inner join cpt_cone_penetrometer_survey \ on cpt_geotechnical_survey.geotechnical_survey_id =cpt_cone_penetrometer_survey.geotechnical_survey_id\ INNER JOIN \ cpt_cone_penetration_test on cpt_cone_penetrometer_survey.cone_penetrometer_survey_id=cpt_cone_penetration_test.cone_penetration_test_id\ INNER JOIN \ cpt_cone_penetration_test_result on cpt_cone_penetration_test_result.cone_penetration_test_id = cpt_cone_penetration_test.cone_penetration_test_id\ INNER JOIN\ cpt_cone_penetrometer on cpt_cone_penetrometer.cone_penetrometer_id = cpt_cone_penetration_test.cone_penetration_test_id\ where cpt_geotechnical_survey.bro_id IN (" for k in range (0,len(cpts)-1): sql=sql + " '%s', "%cpts['bro_id'].iloc[k] sql=sql+"'%s' )"%cpts['bro_id'].iloc[-1] sql=sql+" ORDER BY 'depth'" # print(sql) cursor.execute(sql) self.results = pd.DataFrame(cursor.fetchall(),columns=('penetration_length', 'cone_resistance', 'depth', 'local_friction', 'friction_ratio', 'inclination_resultant', 'pore_pressure_u1', 'pore_pressure_u2', 'pore_pressure_u3', 'x_or_lon', 'y_or_lat', 'offset', 'bro_id', 'cone_surface_quotient' )) with open('cpts.pkl', 'wb') as f: pickle.dump(self.results, f) # plt.scatter(cpts['x_or_lon'].to_numpy(dtype='float'),cpts['y_or_lat'].to_numpy(dtype='float')) self.results.to_csv('cpts.csv') def plot_cpts_one_by_one(self): int2=VtkClass() sel_loc=self.resutls['bro_id'].unique() for i in range(0,len(sel_loc)): data=self.resutls.loc[self.resutls['bro_id']==sel_loc[i]] center=np.asarray([self.resutls['x_or_lon'].iloc[0],self.resutls['y_or_lat'].iloc[0]]).astype("float") elev=self.resutls['offset'].iloc[0] data=data[['depth','local_friction','friction_ratio','pore_pressure_u2']].values data=data[data[:, 0].argsort()] data[np.isnan(data)]=0 # data[:,0]=data[:,0]/1000 int2.make_borehole_as_cube(filename='results\\vtk\\'+sel_loc[i]+'.vtk', data=data, center=center, radius=5, elev=elev) def plot_cpts(self,radius=20): int2=VtkClass() # remove nan cpts = self.results[self.results['depth'].notna()] un=cpts['bro_id'].unique() for i in range(0,len(un)): data=cpts.loc[cpts['bro_id']==un[i]] data=data.sort_values(by=['depth']) bot=data['depth'].values top=np.r_[0,bot[1:]] data=data.fillna(0) data=np.c_[top,bot,np.float32(data['x_or_lon'].values),np.float32(data['y_or_lat'].values),data['offset'].values,data['cone_resistance'].values,data['friction_ratio'],data['local_friction'].values,data['pore_pressure_u2'].values] if i==0: out=data else: out=np.r_[out,data] int2.make_borehole_as_cube_multi('cpts.vtk',out,radius,label=['cone_resistance','local_friction','friction_ratio','pore_pressure']) def classify_cpts(self,radius=20): cpts = self.result() un=cpts['bro_id'].unique() int3=cpt_module()