diff --git a/C_classifier.py b/C_classifier.py new file mode 100644 index 0000000..6c8e01a --- /dev/null +++ b/C_classifier.py @@ -0,0 +1,322 @@ +#------------------------------------------------------------------------------- +# +# Signed Distances Function Interpolator +# ************************************** +# +# This SGeMS plugin interpolates, simulates and classify the signed distance +# function calculated, and creates realizations of an geologic model. +# +# +# AUTHORS: Flávio A. N. Amarante and Roberto Mentzingen Rolo +# +#------------------------------------------------------------------------------- + +#!/bin/python +import sgems +import math +import numpy as np +import random +import copy +from scipy.stats import norm +import sys + +#Creates a randon path given the size of the grid +def random_path(prop): + nodes_not_nan = [] + for i in range(len(prop)): + if not math.isnan(prop[i]): + nodes_not_nan.append(i) + random.shuffle(nodes_not_nan) + return nodes_not_nan + +#Calculates the proportions of variables on a list +def proportion(var, RT): + rock_types =[] + target_prop = [] + for k in range(len(RT)): + target_prop.append(0) + rock_types.append(int(RT[k].split('RT_')[-1])) + rock_types.sort() + var_not_nan = [] + for i in var: + if not math.isnan(i): + var_not_nan.append(i) + for i in range(len(rock_types)): + target_prop[i] = float(var.count(rock_types[i]))/len(var_not_nan) + return target_prop + +def create_variable(grid, name, list): + lst_props_grid = sgems.get_property_list(grid) + prop_final_data_name = name + + if (prop_final_data_name in lst_props_grid): + flag = 0 + i = 1 + while (flag == 0): + test_name = prop_final_data_name + '-' + str(i) + if (test_name not in lst_props_grid): + flag = 1 + prop_final_data_name = test_name + i = i + 1 + + sgems.set_property(grid, prop_final_data_name, list) + +def classifier(grid, sim_val, interpol_val, C, var_n): + #sim_val = sgems.get_property(grid, sim_val) + #interpol_val = sgems.get_property(grid, interpol_val) + + trans_sim = sim_val + for n,i in enumerate(sim_val): + if not math.isnan(i): + trans_sim[n] = (2*C*norm.cdf(i)-C) + + class_blocks = trans_sim + for n,i in enumerate(sim_val): + if not math.isnan(i): + if i > interpol_val[n]: + class_blocks[n] = 1 + if i < interpol_val[n]: + class_blocks[n] = 0 + for n,i in enumerate(interpol_val): + if i < -C: + class_blocks[n] = 1 + if i > C: + class_blocks[n] = 0 + + FN = "Class_C_("+str(var_n)+")_C_"+ str(C) + create_variable(grid, FN, class_blocks) + +#Transform i,j,k in n +def ijk_in_n(grid, i, j, k): + dims = sgems.get_dims(grid) + n = k*dims[0]*dims[1]+j*dims[0]+i + return n + +#Crestes a list with indices of the neighbors valid blocks +def neighb(grid, indice): + ijk = sgems.get_ijk(grid, indice) + neighborhood = [] + for i in range(ijk[0]-1,ijk[0]+2): + for j in range(ijk[1]-1,ijk[1]+2): + for k in range(ijk[2]-1,ijk[2]+2): + ijk_blk = [i,j,k] + neighborhood.append(ijk_blk) + dims = sgems.get_dims(grid) + neighborhood_cp = copy.copy(neighborhood) + for i in neighborhood_cp: + if dims[2] == 1: + if i[0] < 0 or i[1] < 0: + neighborhood.remove(i) + elif i[0] > (dims[0] - 1) or i[1] > (dims[1] - 1): + neighborhood.remove(i) + elif i[2] != 0: + neighborhood.remove(i) + elif i == sgems.get_ijk(grid, indice): + neighborhood.remove(i) + else: + if i[0] < 0 or i[1] < 0 or i[2] < 0: + neighborhood.remove(i) + elif i[0] > (dims[0] - 1) or i[1] > (dims[1] - 1) or i[2] > (dims[2] - 1): + neighborhood.remove(i) + elif i == sgems.get_ijk(grid, indice): + neighborhood.remove(i) + neighborhood_n = [] + for i in neighborhood: + neighborhood_n.append(ijk_in_n(grid,i[0],i[1],i[2])) + return neighborhood_n + + +# Shows every parameter of the plugin in the command pannel +def read_params(a, j=''): + for i in a: + if (type(a[i]) != type({'a': 1})): + print j + "['" + str(i) + "']=" + str(a[i]) + else: + read_params(a[i], j + "['" + str(i) + "']") + +class C_classifier: + def __init__(self): + pass + + def initialize(self, params): + self.params = params + return True + + def execute(self): + + #Execute the funtion read_params + #read_params(self.params) + #print self.params''' + + #temp = sys.stdout + #sys.stdout = sys.stderr + #sys.stderr = temp + + # Get the grid and rock type propery + grid = self.params['propertyselectornoregion']['grid'] + prop = self.params['propertyselectornoregion']['property'] + + # Get the X, Y and Z coordinates and RT property + X = sgems.get_property(grid, '_X_') + Y = sgems.get_property(grid, '_Y_') + Z = sgems.get_property(grid, '_Z_') + RT_data = sgems.get_property(grid, prop) + + # Getting properties + grid_krig = self.params['gridinter']['value'] + grid_var = self.params['gridselectorbasic']['value'] + props = (self.params['orderedpropertyselector']['value']).split(';') + n_var = int(self.params['indicator_regionalization_input']['number_of_indicator_group']) + n_prop = int(self.params['orderedpropertyselector']['count']) + min_cond = self.params['spinBox_2']['value'] + max_cond = self.params['spinBox']['value'] + C_val = float(self.params['c_par']['value']) + + # Error messages + if len(grid_var) == 0 or len(grid_krig) == 0: + print 'Select the variables' + return False + + if n_var != n_prop: + print 'Number of variables and number of variograms models are diferent.' + return False + + # Creating an empty list to store the interpolated distances + SG_OK_list = [] + + # Loop in every variable + for i in xrange(0, n_var): + + #....... Interpolation .........# + + #Getting variables Interpolation + prop_HD = props[i] + prop_name = "Interpolated_" + str(prop_HD) + prop_name_var = "Interpolated_" + str(prop_HD) + ' krig_var' + var_str = '' + indicator_group = "Indicator_group_" + str(i + 1) + elipsoide = self.params['ellipsoidinput']['value'] + n_struct = int(self.params['indicator_regionalization_input'][indicator_group]['Covariance_input']['structures_count']) + + # Error message + if n_struct == 0: + print 'Variogram have no structures' + #continue + return False + + # Loop in every variogram structure - Interpolation + for j in xrange(0, n_struct): + + # Getting variogram parameters + Structure = "Structure_" + str(j + 1) + cov_type = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['type'] + cont = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['contribution'] + + if cov_type == 'Nugget Covariance': + # Writing variogram parameters on a variable in nugget effect case + var_str = var_str + '<{} type="{}"> '.format(Structure, 'Covariance', cont, cov_type, Structure) + print var_str + else: + range1 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range1'] + range2 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range2'] + range3 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range3'] + rake = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['rake'] + dip = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['dip'] + azimuth = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['azimuth'] + + # Writing variogram parameters on a variable in other cases + var_str = var_str + '<{} type="{}"> '.format(Structure, 'Covariance', cont, cov_type, range1, range2, range3, azimuth, dip, rake, Structure) + + # Calling ordinary kriging for each variable, using the variograms parameters above + sgems.execute('RunGeostatAlgorithm kriging::/GeostatParamUtils/XML:: {} '.format(n_struct, var_str, grid_krig, prop_name, grid_var, prop_HD, min_cond, max_cond, elipsoide)) + + SG_OK_list = sgems.get_property(grid_krig, prop_name) + + # Uncertainty Zone + zone_name = "UZ_C_"+str(C_val) + zone = np.array(SG_OK_list) + m1 = zone < C_val + m2 = zone > -C_val + mask = np.logical_and(m1,m2) + uzz = mask.astype('int').tolist() + sgems.set_region(grid_krig,zone_name,uzz) + + #......... Simulation .........# + + #Creating an empty list to store the simulated distances + SG_sim_list = [] + + # Getting variables Simulation + grid_sim_region = zone_name + prop_name_sim = "Sim_C_"+ str(C_val) + var_str_s = '' + nb_realz = self.params['spinBox_3']['value'] + max_cond_hard_data = self.params['Max_Conditioning_Data']['value'] + max_cond_prev_sim_data = self.params['Max_Conditioning_Simul_Data']['value'] + number_pro = self.params['spinBox_6']['value'] + sim_group = "Indicator_group_" + str(i + 1) + elipsoide_sim = self.params['Search_Ellipsoid_Sim']['value'] + n_struct_s = int(self.params['sim_regionalization_input'][sim_group]['Covariance_input']['structures_count']) + HD = "" + p_HD = "" + + #Error message + if n_struct_s == 0: + print 'Variogram have no structures' + return False + + # Loop in every variogram structure + for j in xrange(0, n_struct): + + # Getting variogram parameters + Structure_s = "Structure_" + str(j + 1) + cov_type_s = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['type'] + cont_s = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['contribution'] + + if cov_type_s == 'Nugget Covariance': + #Writing variogram parameters on a variable in nugget effect case + var_str_s = var_str + '<{} type="{}"> '.format(Structure, 'Covariance', cont, cov_type, Structure) + + else: + srange1 = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['ranges']['range1'] + srange2 = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['ranges']['range2'] + srange3 = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['ranges']['range3'] + srake = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['angles']['rake'] + sdip = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['angles']['dip'] + sazimuth = self.params['sim_regionalization_input'][sim_group]['Covariance_input'][Structure_s]['Two_point_model']['angles']['azimuth'] + + # Writing variogram parameters on a variable in other cases + var_str_s = var_str + '<{} type="{}"> '.format(Structure_s, 'Covariance', cont_s, cov_type_s, srange1, srange2, srange3, sazimuth, sdip, srake, Structure_s) + + # Calling sgsim, using the variograms parameters above + sgems.execute('RunGeostatAlgorithm sgsim::/GeostatParamUtils/XML:: {} '.format(number_pro, grid_krig, grid_sim_region, prop_name_sim, nb_realz, max_cond_hard_data, max_cond_prev_sim_data, elipsoide_sim, HD, p_HD, grid_var, prop_HD, n_struct_s, var_str_s)) + + Simulations = sgems.get_properties_in_group(grid_krig,prop_name_sim) + + for prop in Simulations: + ss = sgems.get_property(grid_krig, prop) + SG_sim_list.append(ss) + + #classifying categories + for i in SG_sim_list: + classifier(grid_krig,i,SG_OK_list,C_val,prop_HD) + + #Deleting interpolation and simulations + if self.params['Keep_in']['value']=='0': + sgems.execute('DeleteObjectProperties {}::{}'.format(grid_krig, prop_name)) + if self.params['keep_sim']['value']=='0': + sgems.execute('DeleteObjectProperties {}::{}'.format(grid_krig, prop_name_sim)) + + return True + + def finalize(self): + + return True + + def name(self): + + return "C_classifier" + +################################################################################ +def get_plugins(): + return ["C_classifier"] diff --git a/interpolator.ui b/C_classifier.ui old mode 100755 new mode 100644 similarity index 50% rename from interpolator.ui rename to C_classifier.ui index cfc4f97..92cfa4e --- a/interpolator.ui +++ b/C_classifier.ui @@ -9,7 +9,7 @@ 0 0 - 475 + 714 964 @@ -20,11 +20,14 @@ DFMod - interpolator + C_classifier + + true + 0 @@ -33,6 +36,18 @@ General + + + + Search strategy + + + + + + + + @@ -52,12 +67,12 @@ - + - + Choose signed distances properties @@ -72,7 +87,21 @@ - + + + + + + C Parameter + + + + + + + + + Conditioning data @@ -118,23 +147,11 @@ - - - - Search strategy - - - - - - - - - Variogram + Variogram interpolation @@ -149,103 +166,207 @@ - + - Options + Simulation - - - - - true - - - Softmax transformation - - - - - - - Gamma - - + + + + + + + Nb of realizations + + + + + + + + 0 + 0 + + + + + 0 + 0 + + + + + 16777215 + 16777215 + + + + 99999999 + + + + - - - - false - - - 0.010000000000000 - - - 999999.000000000000000 - - - 175.000000000000000 - - - - - - - false + + + + + 0 + 0 + - - Servo system + + Search Ellipsoid + + + + + + + Max conditioning hard data + + + false + + + + + + + 0 + + + 20000 + + + 12 + + + + + + + Max conditioning previously simulated data + + + false + + + + + + + 0 + + + 20000 + + + 12 + + + + + + + + + - - + + - Lambda + Numer of processors - - - - false - - - 0.990000000000000 - + + - 0.500000000000000 + 0 - - - - Select target proportion - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - + + + label + spinBox_6 + sequential_simulation_parameters_widget + GroupBox9 + + + + Variogram Simulation + + + + + 10 + 30 + 674 + 786 + + + + + + + 10 + 10 + 674 + 15 + + + + Enter simulation variograms + + + + + + Options + + + + true + + + + 10 + 10 + 681 + 21 + + + + Keep Interpolation + + + + + true + + + + 10 + 40 + 681 + 21 + + + + Keep Simulation + + - + Qt::Vertical @@ -274,11 +395,15 @@ OrderedPropertySelector GsTLGroupBox
qtplugins/selectors.h
+ 1 - PropertySelector + GridSelector QWidget
qtplugins/selectors.h
+ + activated(QString) +
Indicator_regionalization_input @@ -294,62 +419,20 @@ QWidget
qtplugins/ellipsoid_input.h
+ + Sequential_simulation_parameters_widget + QWidget +
qtplugins/sequential_simulation_parameters_widget.h
+
GsTLGroupBox QGroupBox
qtplugins/selectors.h
+ 1
- - softmax_check - toggled(bool) - Gamma - setEnabled(bool) - - - 178 - 60 - - - 313 - 87 - - - - - softmax_check - toggled(bool) - servo_check - setEnabled(bool) - - - 97 - 58 - - - 113 - 121 - - - - - servo_check - toggled(bool) - Lambda - setEnabled(bool) - - - 227 - 117 - - - 319 - 147 - - - gridselectorbasic activated(QString) @@ -357,12 +440,12 @@ show_properties(QString) - 278 - 134 + 301 + 278 - 269 - 196 + 292 + 464 diff --git a/README.md b/README.md deleted file mode 100755 index b924800..0000000 --- a/README.md +++ /dev/null @@ -1,5 +0,0 @@ -# dfmod - -This is an implicit geologic modeling SGeMS plugin based on signed distances functions. UI and PY files are included, just paste them in their folders, or run install.sh bash script if you're on linux, then open SGeMS. - -Created by Roberto Mentzingen Rolo diff --git a/interpolator.py b/interpolator.py deleted file mode 100755 index 21226bd..0000000 --- a/interpolator.py +++ /dev/null @@ -1,382 +0,0 @@ -#------------------------------------------------------------------------------- -# -# Signed Distances Function Interpolator -# ************************************** -# -# This SGeMS plugin interpolates (OK) the signed distance function calculated -# for each data and rock type, and creates a geologic model based on the minimum -# estimated distance. -# -# AUTHOR: Roberto Mentzingen Rolo -# -#------------------------------------------------------------------------------- - -#!/bin/python -import sgems -import math -import numpy as np -import random -import copy - -#Creates a randon path given the size of the grid -def random_path(prop): - nodes_not_nan = [] - for i in range(len(prop)): - if not math.isnan(prop[i]): - nodes_not_nan.append(i) - random.shuffle(nodes_not_nan) - return nodes_not_nan - -#Calculates the proportions of variables on a list -def proportion(var, RT): - rock_types =[] - target_prop = [] - for k in range(len(RT)): - target_prop.append(0) - rock_types.append(int(RT[k].split('RT_')[-1])) - rock_types.sort() - var_not_nan = [] - for i in var: - if not math.isnan(i): - var_not_nan.append(i) - for i in range(len(rock_types)): - target_prop[i] = float(var.count(rock_types[i]))/len(var_not_nan) - return target_prop - -#Transform i,j,k in n -def ijk_in_n(grid, i, j, k): - dims = sgems.get_dims(grid) - n = k*dims[0]*dims[1]+j*dims[0]+i - return n - -#Crestes a list with indices of the neighbors valid blocks -def neighb(grid, indice): - ijk = sgems.get_ijk(grid, indice) - neighborhood = [] - for i in range(ijk[0]-1,ijk[0]+2): - for j in range(ijk[1]-1,ijk[1]+2): - for k in range(ijk[2]-1,ijk[2]+2): - ijk_blk = [i,j,k] - neighborhood.append(ijk_blk) - dims = sgems.get_dims(grid) - neighborhood_cp = copy.copy(neighborhood) - for i in neighborhood_cp: - if dims[2] == 1: - if i[0] < 0 or i[1] < 0: - neighborhood.remove(i) - elif i[0] > (dims[0] - 1) or i[1] > (dims[1] - 1): - neighborhood.remove(i) - elif i[2] != 0: - neighborhood.remove(i) - elif i == sgems.get_ijk(grid, indice): - neighborhood.remove(i) - else: - if i[0] < 0 or i[1] < 0 or i[2] < 0: - neighborhood.remove(i) - elif i[0] > (dims[0] - 1) or i[1] > (dims[1] - 1) or i[2] > (dims[2] - 1): - neighborhood.remove(i) - elif i == sgems.get_ijk(grid, indice): - neighborhood.remove(i) - neighborhood_n = [] - for i in neighborhood: - neighborhood_n.append(ijk_in_n(grid,i[0],i[1],i[2])) - return neighborhood_n - - -# Shows every parameter of the plugin in the command pannel -def read_params(a, j=''): - for i in a: - if (type(a[i]) != type({'a': 1})): - print j + "['" + str(i) + "']=" + str(a[i]) - else: - read_params(a[i], j + "['" + str(i) + "']") - -class interpolator: - def __init__(self): - pass - - def initialize(self, params): - self.params = params - return True - - def execute(self): - - '''# Execute the funtion read_params - read_params(self.params) - print self.params''' - - #Get the grid and rock type propery - grid = self.params['propertyselectornoregion']['grid'] - prop = self.params['propertyselectornoregion']['property'] - - #Get the X, Y and Z coordinates and RT property - X = sgems.get_property(grid, '_X_') - Y = sgems.get_property(grid, '_Y_') - Z = sgems.get_property(grid, '_Z_') - RT_data = sgems.get_property(grid, prop) - - # Getting properties - grid_krig = self.params['gridselectorbasic_2']['value'] - grid_var = self.params['gridselectorbasic']['value'] - props = (self.params['orderedpropertyselector']['value']).split(';') - n_var = int(self.params['indicator_regionalization_input']['number_of_indicator_group']) - n_prop = int(self.params['orderedpropertyselector']['count']) - min_cond = self.params['spinBox_2']['value'] - max_cond = self.params['spinBox']['value'] - - # Error messages - if len(grid_var) == 0 or len(grid_krig) == 0: - print 'Select the variables' - return False - - if n_var != n_prop: - print 'Number of variables and number of variograms models are diferent.' - return False - - #Creating an empty list to store the interpolated distances - SG_OK_list = [] - - # Loop in every variable - for i in xrange(0, n_var): - - # Getting variables - prop_HD = props[i] - prop_name = "Interpolated_" + str(prop_HD) - prop_name_var = "Interpolated_" + str(prop_HD) + ' krig_var' - var_str = '' - indicator_group = "Indicator_group_" + str(i + 1) - elipsoide = self.params['ellipsoidinput']['value'] - n_struct = int(self.params['indicator_regionalization_input'][indicator_group]['Covariance_input']['structures_count']) - - # Error message - if n_struct == 0: - print 'Variogram have no structures' - return False - - # Loop in every variogram structure - for j in xrange(0, n_struct): - # Getting variogram parameters - Structure = "Structure_" + str(j + 1) - - cov_type = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['type'] - - cont = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['contribution'] - - if cov_type == 'Nugget Covariance': - #Writing variogram parameters on a variable in nugget effect case - var_str = var_str + '<{} type="{}"> '.format(Structure, 'Covariance', cont, cov_type, Structure) - - else: - range1 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range1'] - range2 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range2'] - range3 = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['ranges']['range3'] - - rake = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['rake'] - dip = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['dip'] - azimuth = self.params['indicator_regionalization_input'][indicator_group]['Covariance_input'][Structure]['Two_point_model']['angles']['azimuth'] - - # Writing variogram parameters on a variable in other cases - var_str = var_str + '<{} type="{}"> '.format(Structure, 'Covariance', cont, cov_type, range1, range2, range3, azimuth, dip, rake, Structure) - - # Calling ordinary kriging for each variable, using the variograms parameters above - sgems.execute('RunGeostatAlgorithm kriging::/GeostatParamUtils/XML:: {} '.format(n_struct, var_str, grid_krig, prop_name, grid_var, prop_HD, min_cond, max_cond, elipsoide)) - - SG_OK_list.append(sgems.get_property(grid_krig, prop_name)) - - #Deleting kriged distances - sgems.execute('DeleteObjectProperties {}::{}'.format(grid_krig, prop_name)) - sgems.execute('DeleteObjectProperties {}::{}'.format(grid_krig, prop_name_var)) - - RT = (self.params['orderedpropertyselector']['value']).split(';') - - #Determinig geomodel based on minimum estimed signed distance function - GeoModel = SG_OK_list[0][:] - - t = 0 - for i in range(len(SG_OK_list[0])): - sgmin = 10e21 - for j in range(len(SG_OK_list)): - if SG_OK_list[j][i] < sgmin: - sgmin = SG_OK_list[j][i] - t = j - if math.isnan(SG_OK_list[j][i]): - GeoModel[i] = float('nan') - else: - GeoModel[i] = (int(RT[t].split('RT_')[-1])) - - #Creating GeoModel property - lst_props_grid=sgems.get_property_list(grid_krig) - prop_final_data_name = 'Geologic_Model' - - if (prop_final_data_name in lst_props_grid): - flag=0 - i=1 - while (flag==0): - test_name=prop_final_data_name+'-'+str(i) - if (test_name not in lst_props_grid): - flag=1 - prop_final_data_name=test_name - i=i+1 - - #Assign conditioning data to grid node - for i in range(len(RT_data)): - if not math.isnan(RT_data[i]): - closest_node = sgems.get_closest_nodeid(grid_krig, X[i],Y[i],Z[i]) - GeoModel[closest_node] = RT_data[i] - - sgems.set_property(grid_krig, prop_final_data_name, GeoModel) - - #Operating softmax transformation - if self.params['softmax_check']['value']=='1': - - gamma =float( self.params['Gamma']['value']) - Prob_list = SG_OK_list[:] - - for i in range(len(SG_OK_list[0])): - soma = 0 - for j in range(len(SG_OK_list)): - soma = soma + math.exp(-SG_OK_list[j][i]/gamma) - for j in range(len(SG_OK_list)): - Prob_list[j][i] = math.exp(-SG_OK_list[j][i]/gamma)/soma - - #Creating probabilities propreties - for k in range(len(Prob_list)): - prop_final_data_name = 'Probability_RT'+str(RT[k].split('RT_')[-1]) - - if (prop_final_data_name in lst_props_grid): - flag=0 - i=1 - while (flag==0): - test_name=prop_final_data_name+'-'+str(i) - if (test_name not in lst_props_grid): - flag=1 - prop_final_data_name=test_name - i=i+1 - - sgems.set_property(grid_krig, prop_final_data_name, Prob_list[k]) - - #Operating servo-system - if self.params['servo_check']['value'] == '1': - var_rt_grid = self.params['targe_prop']['grid'] - var_rt_st = self.params['targe_prop']['property'] - var_rt_region = self.params['targe_prop']['region'] - if len(var_rt_grid) == 0 or len(var_rt_st) == 0: - print 'Select the target proportion property' - return False - - #Getting variables - var_rt = sgems.get_property(var_rt_grid, var_rt_st) - - #Getting parameters - lambda1 = float(self.params['Lambda']['value']) - mi = lambda1/(1-lambda1) - - #Checking if a region exist - if len(var_rt_region) == 0: - #Variable without a region - var_region = var_rt - - else: - region_rt = sgems.get_region(var_rt_grid, var_rt_region) - #Geting the variable inside the region - var_region = [] - for i in range(len(var_rt)): - if region_rt[i] == 1: - var_region.append(var_rt[i]) - - #Getting the target proportion - target_prop = proportion(var_region, RT) - - #Getting the random path - ran_path = random_path(Prob_list[0]) - - #Removing the blocks outside the region from randon path - if len(var_rt_region) != 0: - for i in range(len(region_rt)): - if region_rt[i] == 0: - ran_path.remove(i) - - #servo system - p = 0 - GeoModel_corrected = GeoModel[:] - - visited_rts = [] - for j in ran_path: - visited_rts.append(GeoModel[j]) - instant_proportions = proportion(visited_rts,RT) - - sgmax = 10e-21 - for i in range(len(Prob_list)): - Prob_list[i][j] = Prob_list[i][j] + (mi * (target_prop[i] - instant_proportions[i])) - if Prob_list[i][j] > sgmax: - sgmax = Prob_list[i][j] - p = i - - GeoModel_corrected[j] = int(RT[p][-1]) - visited_rts[-1] = int(RT[p].split('RT_')[-1]) - - #Correcting servo servo-system by the biggest proportion on a neighborhood - GeoModel_corrected_servo_prop = GeoModel_corrected[:] - ran_path_servo_correction = random_path(GeoModel_corrected_servo_prop) - for i in ran_path_servo_correction: - vizinhanca = neighb(grid_krig,i) - - blk_geo_model_corrected_servo = [] - for j in vizinhanca: - blk_geo_model_corrected_servo.append(GeoModel_corrected_servo_prop[j]) - - proportions_servo = proportion(blk_geo_model_corrected_servo, RT) - indice_max_prop = proportions_servo.index(max(proportions_servo)) - - GeoModel_corrected_servo_prop[i] = int(RT[indice_max_prop].split('RT_')[-1]) - - #Creating Geologic_Model_Servo_System property - prop_final_data_name = 'Geologic_Model_Servo_System' - - if (prop_final_data_name in lst_props_grid): - flag=0 - i=1 - while (flag==0): - test_name=prop_final_data_name+'-'+str(i) - if (test_name not in lst_props_grid): - flag=1 - prop_final_data_name=test_name - i=i+1 - - #Creating Geologic_Model_Corrected property - prop_final_data_name1 = 'Geologic_Model_Corrected' - - if (prop_final_data_name1 in lst_props_grid): - flag=0 - i=1 - while (flag==0): - test_name1=prop_final_data_name1+'-'+str(i) - if (test_name1 not in lst_props_grid): - flag=1 - prop_final_data_name1=test_name1 - i=i+1 - - #Assign conditioning data to grid node - for i in range(len(RT_data)): - if not math.isnan(RT_data[i]): - closest_node = sgems.get_closest_nodeid(grid_krig, X[i],Y[i],Z[i]) - GeoModel_corrected[closest_node] = RT_data[i] - GeoModel_corrected_servo_prop[closest_node] = RT_data[i] - - #Setting properties - sgems.set_property(grid_krig, prop_final_data_name, GeoModel_corrected) - sgems.set_property(grid_krig, prop_final_data_name1, GeoModel_corrected_servo_prop) - - return True - - def finalize(self): - - return True - - def name(self): - - return "interpolator" - -################################################################################ -def get_plugins(): - return ["interpolator"] diff --git a/signed_distances.ui b/signed_distances.ui deleted file mode 100755 index 69ed101..0000000 --- a/signed_distances.ui +++ /dev/null @@ -1,83 +0,0 @@ - - - Form - - - true - - - - 0 - 0 - 439 - 536 - - - - Form - - - DFMod - - - signed_distances - - - - - - - - Choose rock property - - - - - - - - - - - - Anisotropic search - - - - - - - - - - - - Qt::Vertical - - - - 20 - 40 - - - - - - - - - - - PropertySelectorNoRegion - QWidget -
qtplugins/selectors.h
-
- - EllipsoidInput - QWidget -
qtplugins/ellipsoid_input.h
-
-
- - -
diff --git a/signed_distances.py b/signed_distances_C.py old mode 100755 new mode 100644 similarity index 90% rename from signed_distances.py rename to signed_distances_C.py index 31f1591..1c32a1e --- a/signed_distances.py +++ b/signed_distances_C.py @@ -1,189 +1,190 @@ -#------------------------------------------------------------------------------- -# -# Signed Distance Function Calculator -# *********************************** -# -# This SGeMS plugin calculates the anisotropic signed distances for each data -# point and each rock type -# -# AUTHOR: Roberto Mentzingen Rolo -# -#------------------------------------------------------------------------------- - -#!/bin/python -import sgems -import math -import numpy as np - -#Calculates the distances -def dist(x1, y1, z1, x2, y2, z2): - return math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2) - -#Defines the rotation and dilatation matrices -def rot(range1, range2, range3, azimuth, dip, rake, vetor): - - if azimuth >= 0 and azimuth <=270: - alpha = math.radians(90-azimuth) - else: - alpha = math.radians(450-azimuth) - beta = -math.radians(dip) - phi = math.radians(rake) - - rot_matrix = np.zeros((3,3)) - - rot_matrix[0,0] = math.cos(beta)*math.cos(alpha) - rot_matrix[0,1] = math.cos(beta)*math.sin(alpha) - rot_matrix[0,2] = -math.sin(beta) - rot_matrix[1,0] = (range1/range2)*(-math.cos(phi)*math.sin(alpha)+math.sin(phi)*math.sin(beta)*math.cos(alpha)) - rot_matrix[1,1] = (range1/range2)*(math.cos(phi)*math.cos(alpha)+math.sin(phi)*math.sin(beta)*math.sin(alpha)) - rot_matrix[1,2] = (range1/range2)*(math.sin(phi)*math.cos(beta)) - rot_matrix[2,0] = (range1/range3)*(math.sin(phi)*math.sin(alpha)+math.cos(phi)*math.sin(beta)*math.cos(alpha)) - rot_matrix[2,1] = (range1/range3)*(-math.sin(phi)*math.cos(alpha)+math.cos(phi)*math.sin(beta)*math.sin(alpha)) - rot_matrix[2,2] = (range1/range3)*(math.cos(phi)*math.cos(beta)) - - vetor = np.array(vetor) - - return np.dot(rot_matrix, vetor) - -#Transform the data with the ratation/dilatation matrices -def anis_search(X, Y, Z, range1, range2, range3, azimuth, dip, rake): - - X_linha = [] - Y_linha = [] - Z_linha = [] - - for i in range(len(X)): - vet = [X[i],Y[i],Z[i]] - - vet_rot = rot(range1, range2, range3, azimuth, dip, rake, vet) - - X_linha.append(vet_rot[0]) - Y_linha.append(vet_rot[1]) - Z_linha.append(vet_rot[2]) - - return X_linha, Y_linha, Z_linha - -#Shows every parameter of the plugin in the command pannel -def read_params(a,j=''): - for i in a: - if (type(a[i])!=type({'a':1})): - print j+"['"+str(i)+"']="+str(a[i]) - else: - read_params(a[i],j+"['"+str(i)+"']") - -class signed_distances: - def __init__(self): - pass - - def initialize(self, params): - self.params = params - return True - - def execute(self): - - '''#Execute the funtion read_params - read_params(self.params) - print self.params''' - - #Get the grid and rock type propery - grid = self.params['propertyselectornoregion']['grid'] - prop = self.params['propertyselectornoregion']['property'] - - #Error message - if len(grid) == 0 or len(prop) == 0: - print 'Select the rocktype property' - return False - - #Get the X, Y and Z coordinates and RT property - X = sgems.get_property(grid, '_X_') - Y = sgems.get_property(grid, '_Y_') - Z = sgems.get_property(grid, '_Z_') - RT = sgems.get_property(grid, prop) - - elipsoide = self.params['ellipsoidinput']['value'] - elipsoide_split = elipsoide.split() - - range1 = float(elipsoide_split[0]) - range2 = float(elipsoide_split[1]) - range3 = float(elipsoide_split[2]) - - azimuth = float(elipsoide_split[3]) - dip = float(elipsoide_split[4]) - rake = float(elipsoide_split[5]) - - X, Y, Z = anis_search(X, Y, Z, range1, range2, range3, azimuth, dip, rake) - - #Creates a list of all rock types - rt_list = [] - for i in RT: - if i not in rt_list and not math.isnan(i): - rt_list.append(i) - - #Sort the rock type list in crescent order - rt_list = [int(x) for x in rt_list] - rt_list.sort() - - #Create a empty distance matrix - dist_matrix = np.zeros(shape = ((len(rt_list)), (len(RT)))) - - #Calculates the signed distances, and append it in the distance matrix - for i in range(len(rt_list)): - rock = rt_list[i] - - for j in range(len(RT)): - - if math.isnan(RT[j]): - dist_matrix[i][j] = float('nan') - - elif RT[j] == rock: - dsmin = 1.0e21 - - for k in range(len(RT)): - - if RT[j] != RT[k] and not math.isnan(RT[k]): - if (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) < dsmin: - dsmin = (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) - - dist_matrix[i][j] = -dsmin - - else: - dsmin = 1.0e21 - - for k in range(len(RT)): - - if RT[k] == rock: - if (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) < dsmin: - dsmin = (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) - - dist_matrix[i][j] = dsmin - - #Creates the signed distances properties - lst_props_grid=sgems.get_property_list(grid) - - for k in range(len(dist_matrix)): - prop_final_data_name = 'Signed_Distances_RT_' + str(rt_list[k]) - - if (prop_final_data_name in lst_props_grid): - flag=0 - i=1 - while (flag==0): - test_name=prop_final_data_name+'-'+str(i) - if (test_name not in lst_props_grid): - flag=1 - prop_final_data_name=test_name - i=i+1 - - list = dist_matrix[k].tolist() - sgems.set_property(grid, prop_final_data_name, list) - - return True - - def finalize(self): - return True - - def name(self): - return "signed_distances" - -################################################################################ -def get_plugins(): - return ["signed_distances"] +#------------------------------------------------------------------------------- +# +# Signed Distance Function Calculator C +# ************************************** +# +# This SGeMS plugin calculates the anisotropic signed distances for each data +# point and each rock type +# +# AUTHOR: Roberto Mentzingen Rolo +# Modified by: Flávio A. N. Amarante +#------------------------------------------------------------------------------- + +#!/bin/python +import sgems +import math +import numpy as np + +#Calculates the distances +def dist(x1, y1, z1, x2, y2, z2): + return math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2) + +#Defines the rotation and dilatation matrices +def rot(range1, range2, range3, azimuth, dip, rake, vetor): + + if azimuth >= 0 and azimuth <=270: + alpha = math.radians(90-azimuth) + else: + alpha = math.radians(450-azimuth) + beta = -math.radians(dip) + phi = math.radians(rake) + + rot_matrix = np.zeros((3,3)) + + rot_matrix[0,0] = math.cos(beta)*math.cos(alpha) + rot_matrix[0,1] = math.cos(beta)*math.sin(alpha) + rot_matrix[0,2] = -math.sin(beta) + rot_matrix[1,0] = (range1/range2)*(-math.cos(phi)*math.sin(alpha)+math.sin(phi)*math.sin(beta)*math.cos(alpha)) + rot_matrix[1,1] = (range1/range2)*(math.cos(phi)*math.cos(alpha)+math.sin(phi)*math.sin(beta)*math.sin(alpha)) + rot_matrix[1,2] = (range1/range2)*(math.sin(phi)*math.cos(beta)) + rot_matrix[2,0] = (range1/range3)*(math.sin(phi)*math.sin(alpha)+math.cos(phi)*math.sin(beta)*math.cos(alpha)) + rot_matrix[2,1] = (range1/range3)*(-math.sin(phi)*math.cos(alpha)+math.cos(phi)*math.sin(beta)*math.sin(alpha)) + rot_matrix[2,2] = (range1/range3)*(math.cos(phi)*math.cos(beta)) + + vetor = np.array(vetor) + + return np.dot(rot_matrix, vetor) + +#Transform the data with the ratation/dilatation matrices +def anis_search(X, Y, Z, range1, range2, range3, azimuth, dip, rake): + + X_linha = [] + Y_linha = [] + Z_linha = [] + + for i in range(len(X)): + vet = [X[i],Y[i],Z[i]] + + vet_rot = rot(range1, range2, range3, azimuth, dip, rake, vet) + + X_linha.append(vet_rot[0]) + Y_linha.append(vet_rot[1]) + Z_linha.append(vet_rot[2]) + + return X_linha, Y_linha, Z_linha + +#Shows every parameter of the plugin in the command pannel +def read_params(a,j=''): + for i in a: + if (type(a[i])!=type({'a':1})): + print j+"['"+str(i)+"']="+str(a[i]) + else: + read_params(a[i],j+"['"+str(i)+"']") + +class signed_distances_C: + def __init__(self): + pass + + def initialize(self, params): + self.params = params + return True + + def execute(self): + + '''#Execute the funtion read_params + read_params(self.params) + print self.params''' + + #Get the grid, rock type propery and C value + grid = self.params['propertyselectornoregion']['grid'] + prop = self.params['propertyselectornoregion']['property'] + C = float(self.params['c_par']['value']) + + #Error message + if len(grid) == 0 or len(prop) == 0: + print 'Select the rocktype property' + return False + + #Get the X, Y and Z coordinates and RT property + X = sgems.get_property(grid, '_X_') + Y = sgems.get_property(grid, '_Y_') + Z = sgems.get_property(grid, '_Z_') + RT = sgems.get_property(grid, prop) + + elipsoide = self.params['ellipsoidinput']['value'] + elipsoide_split = elipsoide.split() + + range1 = float(elipsoide_split[0]) + range2 = float(elipsoide_split[1]) + range3 = float(elipsoide_split[2]) + + azimuth = float(elipsoide_split[3]) + dip = float(elipsoide_split[4]) + rake = float(elipsoide_split[5]) + + X, Y, Z = anis_search(X, Y, Z, range1, range2, range3, azimuth, dip, rake) + + #Creates a list of all rock types + rt_list = [] + for i in RT: + if i not in rt_list and not math.isnan(i): + rt_list.append(i) + + #Sort the rock type list in crescent order + rt_list = [int(x) for x in rt_list] + rt_list.sort() + + #Create a empty distance matrix + dist_matrix = np.zeros(shape = ((len(rt_list)), (len(RT)))) + + #Calculates the signed distances, and append it in the distance matrix + for i in range(len(rt_list)): + rock = rt_list[i] + + for j in range(len(RT)): + + if math.isnan(RT[j]): + dist_matrix[i][j] = float('nan') + + elif RT[j] == rock: + dsmin = 1.0e21 + + for k in range(len(RT)): + + if RT[j] != RT[k] and not math.isnan(RT[k]): + if (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) < dsmin: + dsmin = (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) + + dist_matrix[i][j] = -dsmin-C + + else: + dsmin = 1.0e21 + + for k in range(len(RT)): + + if RT[k] == rock: + if (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) < dsmin: + dsmin = (dist(X[j], Y[j], Z[j], X[k], Y[k], Z[k])) + + dist_matrix[i][j] = dsmin+C + + #Creates the signed distances properties + lst_props_grid=sgems.get_property_list(grid) + + for k in range(len(dist_matrix)): + prop_final_data_name = 'Signed_Distances_C_'+ str(C)+'_RT_'+ str(rt_list[k]) + + if (prop_final_data_name in lst_props_grid): + flag=0 + i=1 + while (flag==0): + test_name=prop_final_data_name+'-'+str(i) + if (test_name not in lst_props_grid): + flag=1 + prop_final_data_name=test_name + i=i+1 + + list = dist_matrix[k].tolist() + sgems.set_property(grid, prop_final_data_name, list) + + return True + + def finalize(self): + return True + + def name(self): + return "signed_distances_C" + +################################################################################ +def get_plugins(): + return ["signed_distances_C"] diff --git a/signed_distances_C.ui b/signed_distances_C.ui new file mode 100644 index 0000000..3f794a3 --- /dev/null +++ b/signed_distances_C.ui @@ -0,0 +1,93 @@ + + + Form + + + true + + + + 0 + 0 + 453 + 536 + + + + Form + + + DFMod + + + signed_distances_C + + + + + + Choose rock property + + + + + + + + + + + + Anisotropic search + + + + + + + + + + + + + + C Parameter + + + + + + + + + + + + Qt::Vertical + + + + 20 + 386 + + + + + + + + + PropertySelectorNoRegion + QWidget +
qtplugins/selectors.h
+
+ + EllipsoidInput + QWidget +
qtplugins/ellipsoid_input.h
+
+
+ + +