From 41fa700e4041ee618cbf0fbcc75c228301b56d1d Mon Sep 17 00:00:00 2001 From: nisharam7 <54867844+nisharam7@users.noreply.github.com> Date: Thu, 2 Apr 2020 16:19:50 -0500 Subject: [PATCH 01/15] new comment --- geostatspy/geostats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 51103ed..bc3c1fa 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -38,7 +38,7 @@ def backtr(df,vcol,vr,vrg,zmin,zmax,ltail,ltpar,utail,utpar): :param utpar: upper tail extrapolation parameter :return: TODO """ - + #comment EPSLON=1.0e-20 nd = len(df); nt = len(vr) # number of data to transform and number of data in table backtr = np.zeros(nd) From 63c52a395272e9f32b9989ad1cdbec343b5f4784 Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Wed, 15 Apr 2020 11:08:53 -0500 Subject: [PATCH 02/15] some setup work for kt3d --- geostatspy/geostats.py | 112 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 105 insertions(+), 7 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 51103ed..0087306 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2196,6 +2196,101 @@ def nscore( return ns, vr, wt_ns +# *** all parameters have been copied over form kb2d with respective z parameter added if necessary +# we can probably delete a couple of these if we don't do block kriging +def kt3d ( + df, + xcol, + ycol, + zcol, + vcol, + tmin, + tmax, + nx, + xmn, + xsiz, + ny, + ymn, + ysiz, + nz, + zmn, + zsiz, + nxdis, + nydis, + nzdis, + ndmin, + ndmax, + radius, + ktype, + skmean, + vario, +): + UNEST = -999. + EPSLON = 1.0e-10 + # VERSION = 2.907 + # first = True + PMX = 9999.0 + # MAXSAM = ndmax + 1 + # MAXDIS = nxdis * nydis + # MAXKD = MAXSAM + 1 + # MAXKRG = MAXKD * MAXKD + + # load the data + # trim values outside tmin and tmax + df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] + nd = len(df_extract) + ndmax = min(ndmax,nd) + x = df_extract[xcol].values + y = df_extract[ycol].values + z = df_extract[zcol].values + vr = df_extract[vcol].values + + # set up tree for nearest neighbor search + dp = list((z[i], y[i], x[i]) for i in range(0,nd)) + data_locs = np.column_stack((z, y, x)) + tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) + + # Summary statistics for the data after trimming + avg = vr.mean() + stdev = vr.std() + ss = stdev**2.0 # variance + vrmin = vr.min() + vrmax = vr.max() + + nst = vario['nst'] # num structures + cc = np.zeros(nst) # variance contribution + aa = np.zeros(nst) # major range + it = np.zeros(nst) # type vario structure + aa = np.zeros(nst) # major range + ang_azi = np.zeros(nst) # azimuth + ang_dip = np.zeros(nst) # dip + anis = np.zeros(nst) # anistropy ratio w/ minor range + anis_v = np.zeros(nst) # anistropy ratio w/ vertical range + + # only works w/ nst == 1 or nst == 2 + c0 = vario['nug'] + it[0] = vario['it1'] + cc[0] = vario['cc1'] + ang_azi[0] = vario['azi1'] + ang_dip[0] = vario['dip1'] + aa[0] = vario['hmax'] + anis[0] = vario['hmed1']/vario['hmax1'] + anis_v[0] = vario['hmin1']/vario['hmax1'] + if nst == 2: + cc[1] = vario['cc2'] + it[1] = vario['it2'] + ang_azi[1] = vario['azi2'] + ang_dip[1] = vario['dip1'] + aa[1] = vario['hmax'] + anis[1] = vario['hmed2']/vario['hmax2'] + anis_v[1] = vario['hmin2']/vario['hmax2'] + + rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX) + cbb = maxcov + + # need to set up numPy cube grid + + # main loop over all points: def kb2d( df, @@ -2258,8 +2353,11 @@ def kb2d( # load the variogram nst = vario['nst'] - cc = np.zeros(nst); aa = np.zeros(nst); it = np.zeros(nst) - ang = np.zeros(nst); anis = np.zeros(nst) + cc = np.zeros(nst); + aa = np.zeros(nst); + it = np.zeros(nst) + ang = np.zeros(nst); + anis = np.zeros(nst) c0 = vario['nug']; cc[0] = vario['cc1']; it[0] = vario['it1']; ang[0] = vario['azi1']; @@ -2267,6 +2365,7 @@ def kb2d( if nst == 2: cc[1] = vario['cc2']; it[1] = vario['it2']; ang[1] = vario['azi2']; aa[1] = vario['hmaj2']; anis[1] = vario['hmin2']/vario['hmaj2']; + # Allocate the needed memory: xdb = np.zeros(MAXDIS) @@ -2299,7 +2398,7 @@ def kb2d( # Summary statistics for the data after trimming avg = vr.mean() stdev = vr.std() - ss = stdev**2.0 + ss = stdev**2.0 # variance vrmin = vr.min() vrmax = vr.max() @@ -2324,6 +2423,7 @@ def kb2d( xdb[i] = xloc ydb[i] = yloc + # Initialize accumulators: cbb = 0.0 rad2 = radius*radius @@ -2353,7 +2453,7 @@ def kb2d( yloc = ymn + (iy-0)*ysiz for ix in range(0,nx): xloc = xmn + (ix-0)*xsiz - current_node = (yloc,xloc) + current_node = (yloc,xloc) # xloc, yloc centroid of cell # Find the nearest samples within each octant: First initialize # the counter arrays: @@ -2429,6 +2529,7 @@ def kb2d( # Establish Right Hand Side Covariance: if ndb <= 1: + # ndb <= 1 -> use single value cb = cova2(xx,yy,xdb[0],ydb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) else: cb = 0.0 @@ -4470,9 +4571,6 @@ def cova3(x1, y1, z1, x2, y2, z2, nst, c0, pmx, cc, aa, it, anis, anis_v, rotmat cova3_ = cova3_ + cov1 return cova3_ - - - def gamv_3D( df, From 70057571015db735a66b06b47b7a4e20dd25da54 Mon Sep 17 00:00:00 2001 From: Dale Kang Date: Mon, 20 Apr 2020 22:58:17 -0500 Subject: [PATCH 03/15] setup loops, allocated memory for variables, setup array for value of neighbors --- geostatspy/geostats.py | 50 +++++++++++++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index daf6cc5..1c2f509 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2235,6 +2235,23 @@ def kt3d ( # MAXKD = MAXSAM + 1 # MAXKRG = MAXKD * MAXKD +# Allocate the needed memory: + xdb = np.zeros(MAXDIS) + ydb = np.zeros(MAXDIS) + zdb = np.zeros(MAXDIS) + xa = np.zeros(MAXSAM) + ya = np.zeros(MAXSAM) + za = np.zeros(MAXSAM) + vra = np.zeros(MAXSAM) + dist = np.zeros(MAXSAM) + nums = np.zeros(MAXSAM) + r = np.zeros(MAXKD) + rr = np.zeros(MAXKD) + s = np.zeros(MAXKD) + a = np.zeros(MAXKRG) + kmap = np.zeros((nx,ny,nz)) + vmap = np.zeros((nx,ny,nz)) + # load the data # trim values outside tmin and tmax df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] @@ -2257,6 +2274,7 @@ def kt3d ( vrmin = vr.min() vrmax = vr.max() + # load the variogram nst = vario['nst'] # num structures cc = np.zeros(nst) # variance contribution aa = np.zeros(nst) # major range @@ -2288,9 +2306,39 @@ def kt3d ( rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX) cbb = maxcov + # TODO Set up discretization points per block # need to set up numPy cube grid - # main loop over all points: + # MAIN LOOP OVER ALL THE BLOCKS IN THE GRID: + for iz in range(0,nz): + zloc = zmn + (iz-0)*zsiz + for iy in range(0,ny): + yloc = ymn + (iy-0)*ysiz + for ix in range(0,nx): + xloc = xmn + (ix-0)*xsiz + current_node = (zloc,yloc,xloc) # centroid + + # Find the nearest samples within each octant: First initializethe counter arrays: + na = -1 # accounting for 0 as first index + dist.fill(1.0e+20) + nums.fill(-1) + dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search + # remove any data outside search radius + na = len(dist) + nums = nums[dist Date: Tue, 21 Apr 2020 14:25:42 -0500 Subject: [PATCH 04/15] setup array holding neighborhood samples and situation of only one sample --- geostatspy/geostats.py | 58 +++++++++++++++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 12 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 1c2f509..389b5ea 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2252,7 +2252,7 @@ def kt3d ( kmap = np.zeros((nx,ny,nz)) vmap = np.zeros((nx,ny,nz)) - # load the data +# load the data # trim values outside tmin and tmax df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] nd = len(df_extract) @@ -2262,19 +2262,19 @@ def kt3d ( z = df_extract[zcol].values vr = df_extract[vcol].values - # set up tree for nearest neighbor search +# set up tree for nearest neighbor search dp = list((z[i], y[i], x[i]) for i in range(0,nd)) data_locs = np.column_stack((z, y, x)) tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) - # Summary statistics for the data after trimming +# Summary statistics for the data after trimming avg = vr.mean() stdev = vr.std() ss = stdev**2.0 # variance vrmin = vr.min() vrmax = vr.max() - # load the variogram +# load the variogram nst = vario['nst'] # num structures cc = np.zeros(nst) # variance contribution aa = np.zeros(nst) # major range @@ -2309,7 +2309,7 @@ def kt3d ( # TODO Set up discretization points per block # need to set up numPy cube grid - # MAIN LOOP OVER ALL THE BLOCKS IN THE GRID: +# MAIN LOOP OVER ALL THE BLOCKS IN THE GRID: for iz in range(0,nz): zloc = zmn + (iz-0)*zsiz for iy in range(0,ny): @@ -2318,7 +2318,7 @@ def kt3d ( xloc = xmn + (ix-0)*xsiz current_node = (zloc,yloc,xloc) # centroid - # Find the nearest samples within each octant: First initializethe counter arrays: +# Find the nearest samples within each octant: First initializethe counter arrays: na = -1 # accounting for 0 as first index dist.fill(1.0e+20) nums.fill(-1) @@ -2329,14 +2329,48 @@ def kt3d ( dist = dist[dist Date: Tue, 21 Apr 2020 16:07:53 -0500 Subject: [PATCH 05/15] filling up matrices/creating numpy cube --- geostatspy/geostats.py | 113 ++++++++++++++++++++++++++++++++++++++--- 1 file changed, 106 insertions(+), 7 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index daf6cc5..02eb4dd 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2230,10 +2230,10 @@ def kt3d ( # VERSION = 2.907 # first = True PMX = 9999.0 - # MAXSAM = ndmax + 1 - # MAXDIS = nxdis * nydis - # MAXKD = MAXSAM + 1 - # MAXKRG = MAXKD * MAXKD + MAXSAM = ndmax + 1 + MAXDIS = nxdis * nydis * nzdis + MAXKD = MAXSAM + 1 + MAXKRG = MAXKD * MAXKD # load the data # trim values outside tmin and tmax @@ -2245,10 +2245,27 @@ def kt3d ( z = df_extract[zcol].values vr = df_extract[vcol].values + # Allocate the needed memory: + xdb = np.zeros(MAXDIS) + ydb = np.zeros(MAXDIS) + zdb = np.zeros(MAXDIS) + xa = np.zeros(MAXSAM) + ya = np.zeros(MAXSAM) + za = np.zeros(MAXSAM) + vra = np.zeros(MAXSAM) + dist = np.zeros(MAXSAM) + nums = np.zeros(MAXSAM) + r = np.zeros(MAXKD) + rr = np.zeros(MAXKD) + s = np.zeros(MAXKD) + a = np.zeros(MAXKRG) + kmap = np.zeros((nx,ny,nz)) + vmap = np.zeros((nx,ny,nz)) + # set up tree for nearest neighbor search dp = list((z[i], y[i], x[i]) for i in range(0,nd)) data_locs = np.column_stack((z, y, x)) - tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) + tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error? # Summary statistics for the data after trimming avg = vr.mean() @@ -2289,9 +2306,91 @@ def kt3d ( cbb = maxcov # need to set up numPy cube grid + cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') # main loop over all points: + nk = 0 + ak = 0.0 + vk = 0.0 + for iz in range(0,nz): + zloc = zmn + (iz-0)*zsiz #am I supposed to add these values in the numpy grid? + for iy in range(0,ny): + yloc = ymn + (iy-0)*ysiz + for ix in range (0,nz): #triple nested for loop that loops through points amd adds to matrices + xloc = xmn + (ix-0)*xsiz + current_node = (zloc, yloc,xloc) # xloc, yloc centroid of cell + + # Find the nearest samples within each octant: First initialize +# the counter arrays: + na = -1 # accounting for 0 as first index + dist.fill(1.0e+20) + nums.fill(-1) + dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search + # remove any data outside search radius + na = len(dist) + nums = nums[dist use single value + cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) def kb2d( df, xcol, @@ -2521,7 +2620,7 @@ def kb2d( for i in range(0,na): # was j - want full matrix iin = iin + 1 a[iin] = cova2(xa[i],ya[i],xa[j],ya[j],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) - if ktype == 1: + if ktype == 1: iin = iin + 1 a[iin] = unbias xx = xa[j] - xloc @@ -2531,7 +2630,7 @@ def kb2d( if ndb <= 1: # ndb <= 1 -> use single value cb = cova2(xx,yy,xdb[0],ydb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) - else: + else: cb = 0.0 for j1 in range(0,ndb): cb = cb + cova2(xx,yy,xdb[j1],ydb[j1],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) From 413bf8fae77686853fab175781d35c791bde5a2e Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Tue, 21 Apr 2020 16:25:18 -0500 Subject: [PATCH 06/15] t --- geostatspy/geostats.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 389b5ea..b561150 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2336,7 +2336,7 @@ def kt3d ( print('UNEST at ' + str(ix) + ',' + str(iy) + ',' + str(iz)) else: -# Put coordinates and values of neighborhood samples into xa,ya,vra: +# Put coordinates and values of neighborhood samples into xa,ya,za,vra: for ia in range(0,na): jj = int(nums[ia]) xa[ia] = x[jj] From 1d54ebbd02e197fb8de68dfd5166e66d12e2d9d3 Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Tue, 21 Apr 2020 16:37:47 -0500 Subject: [PATCH 07/15] small fixes --- geostatspy/geostats.py | 81 +++++++++++++++++++++--------------------- 1 file changed, 41 insertions(+), 40 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index a044f5d..8986f90 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2278,7 +2278,6 @@ def kt3d ( kmap = np.zeros((nx,ny,nz)) vmap = np.zeros((nx,ny,nz)) - # set up tree for nearest neighbor search # set up tree for nearest neighbor search dp = list((z[i], y[i], x[i]) for i in range(0,nd)) data_locs = np.column_stack((z, y, x)) @@ -2323,7 +2322,8 @@ def kt3d ( rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX) cbb = maxcov - # TODO Set up discretization points per block + # TODO Set up discretization points per block (bock kriging) + # need to set up numPy cube grid cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') # main loop over all points: @@ -2331,12 +2331,13 @@ def kt3d ( ak = 0.0 vk = 0.0 for iz in range(0,nz): - zloc = zmn + (iz-0)*zsiz #am I supposed to add these values in the numpy grid? + zloc = zmn + (iz-0)*zsiz for iy in range(0,ny): yloc = ymn + (iy-0)*ysiz - for ix in range (0,nz): #triple nested for loop that loops through points amd adds to matrices + for ix in range (0,nz): #triple nested for loop that loops through points and adds to matrices xloc = xmn + (ix-0)*xsiz - current_node = (zloc, yloc,xloc) # xloc, yloc centroid of cell + current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell + # Find the nearest samples within each octant: First initialize # the counter arrays: @@ -2350,14 +2351,14 @@ def kt3d ( dist = dist[dist use single value - cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[0],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) def kb2d( df, xcol, From 43a20201b1313c83903557c6423df538859ceb1f Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Tue, 28 Apr 2020 16:31:13 -0500 Subject: [PATCH 08/15] right hand covariance --- geostatspy/geostats.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 8986f90..cf9e967 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2339,8 +2339,7 @@ def kt3d ( current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell - # Find the nearest samples within each octant: First initialize -# the counter arrays: + # Find the nearest samples within each octant: First initializ the counter arrays: na = -1 # accounting for 0 as first index dist.fill(1.0e+20) nums.fill(-1) @@ -2357,7 +2356,6 @@ def kt3d ( estv = UNEST print('UNEST at ' + str(ix) + ',' + str(iy) + ','+str(iz)) else: - # Put coordinates and values of neighborhood samples into xa,ya,za,vra: for ia in range(0,na): jj = int(nums[ia]) @@ -2377,8 +2375,7 @@ def kt3d ( # Establish Right Hand Side Covariance: if ndb <= 1: cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) - else: - + else: cb = 0.0 for i in range(0,ndb): cb = cb + cova3(xx,yy,zz,xdb[i],ydb[i],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) @@ -2410,6 +2407,7 @@ def kt3d ( for j in range(0,na): for i in range(0,na): #loops through and adds covariances into matrix a[j][i] = cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) + r[j] = cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) # TODO Set up kriging matrices: From f063ead4212c81fbe21c31e79f9f617ed8990520 Mon Sep 17 00:00:00 2001 From: Saahithi Joopelli Date: Wed, 29 Apr 2020 19:28:43 -0500 Subject: [PATCH 09/15] Update geostats.py added in the estimation variance and filling out kmap and vmap --- geostatspy/geostats.py | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index cf9e967..6a0f509 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2409,7 +2409,32 @@ def kt3d ( a[j][i] = cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) r[j] = cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) +# Write a warning if the matrix is singular: + if ising != 0: + print('WARNING KT3D: singular matrix') + print(' for block' + str(ix) + ',' + str(iy) + ',' + str(iz)) + est = UNEST + estv = UNEST + else: +# Compute the estimate and the kriging variance: + est = 0.0 + estv = cbb + sumw = 0.0 + if ktype == 1: + estv = estv - (s[na])*unbias + for i in range(0,na): + sumw = sumw + s[i] + est = est + s[i]*vra[i] + estv = estv - s[i]*rr[i] + if ktype == 0: + est = est + (1.0-sumw)*skmean + kmap[nz-iz-1, ny-iy-1,ix] = est + vmap[nz-iz-1, ny-iy-1,ix] = estv + if est > UNEST: + nk = nk + 1 + ak = ak + est + vk = vk + est*est # TODO Set up kriging matrices: return kmap, vmap From d999def9f8f4de164d47bfd9ec412b0129892c68 Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Tue, 5 May 2020 15:50:37 -0500 Subject: [PATCH 10/15] small fixes --- geostatspy/geostats.py | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 6a0f509..9c680ca 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2321,11 +2321,13 @@ def kt3d ( rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX) cbb = maxcov + unbias = maxcov - # TODO Set up discretization points per block (bock kriging) + # limited to simple kriging for now + ktype = 0 + + # TODO Set up discretization points per block (block kriging) - # need to set up numPy cube grid - cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') # main loop over all points: nk = 0 ak = 0.0 @@ -2409,33 +2411,35 @@ def kt3d ( a[j][i] = cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) r[j] = cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) -# Write a warning if the matrix is singular: + s = ksol_numpy(neq,a,r) + ising = 0 + + # Write a warning if the matrix is singular: if ising != 0: print('WARNING KT3D: singular matrix') - print(' for block' + str(ix) + ',' + str(iy) + ',' + str(iz)) est = UNEST estv = UNEST else: -# Compute the estimate and the kriging variance: + # Compute the estimate and the kriging variance: est = 0.0 estv = cbb sumw = 0.0 - if ktype == 1: - estv = estv - (s[na])*unbias + # if ktype == 1: + # estv = estv - (s[na])*unbias for i in range(0,na): sumw = sumw + s[i] est = est + s[i]*vra[i] estv = estv - s[i]*rr[i] if ktype == 0: est = est + (1.0-sumw)*skmean + # add estimation variance and estimated value to map kmap[nz-iz-1, ny-iy-1,ix] = est vmap[nz-iz-1, ny-iy-1,ix] = estv if est > UNEST: nk = nk + 1 ak = ak + est vk = vk + est*est -# TODO Set up kriging matrices: return kmap, vmap From a7f1478d9f168e6702557988a8a384b40b56e62c Mon Sep 17 00:00:00 2001 From: thingimajig13 <54867822+thingimajig13@users.noreply.github.com> Date: Tue, 5 May 2020 16:06:32 -0500 Subject: [PATCH 11/15] Add files via upload --- geostatspy/testing.ipynb | 471 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 471 insertions(+) create mode 100644 geostatspy/testing.ipynb diff --git a/geostatspy/testing.ipynb b/geostatspy/testing.ipynb new file mode 100644 index 0000000..fc3e27a --- /dev/null +++ b/geostatspy/testing.ipynb @@ -0,0 +1,471 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 88, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "testing the kt3d function using some data from Geostats Guy\n" + ] + } + ], + "source": [ + "print(\"testing the kt3d function using some data from Geostats Guy\")" + ] + }, + { + "cell_type": "code", + "execution_count": 98, + "metadata": {}, + "outputs": [ + { + "ename": "IndentationError", + "evalue": "unindent does not match any outer indentation level (geostats.py, line 364)", + "output_type": "error", + "traceback": [ + "Traceback \u001b[1;36m(most recent call last)\u001b[0m:\n", + " File \u001b[0;32m\"C:\\Users\\meena\\Anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py\"\u001b[0m, line \u001b[0;32m3325\u001b[0m, in \u001b[0;35mrun_code\u001b[0m\n exec(code_obj, self.user_global_ns, self.user_ns)\n", + "\u001b[1;36m File \u001b[1;32m\"\"\u001b[1;36m, line \u001b[1;32m3\u001b[1;36m, in \u001b[1;35m\u001b[1;36m\u001b[0m\n\u001b[1;33m import geostatspy.geostats as geostats\u001b[0m\n", + "\u001b[1;36m File \u001b[1;32m\"C:\\Users\\meena\\AppData\\Roaming\\Python\\Python37\\site-packages\\geostatspy\\geostats.py\"\u001b[1;36m, line \u001b[1;32m364\u001b[0m\n\u001b[1;33m EPSLON=1.0e-20\u001b[0m\n\u001b[1;37m ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m unindent does not match any outer indentation level\n" + ] + } + ], + "source": [ + "import geostatspy.GSLIB as GSLIB \n", + "# GSLIB utilies, visualization and wrapper\n", + "import geostatspy.geostats as geostats\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import scipy.spatial as sp # for fast nearest neighbor search" + ] + }, + { + "cell_type": "code", + "execution_count": 99, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "\n", + "\n", + "# *** all parameters have been copied over form kb2d with respective z parameter added if necessary\n", + "# we can probably delete a couple of these if we don't do block kriging\n", + "def kt3d (\n", + " df,\n", + " xcol,\n", + " ycol,\n", + " zcol,\n", + " vcol,\n", + " tmin,\n", + " tmax,\n", + " nx,\n", + " xmn,\n", + " xsiz,\n", + " ny,\n", + " ymn,\n", + " ysiz,\n", + " nz,\n", + " zmn,\n", + " zsiz,\n", + " nxdis,\n", + " nydis,\n", + " nzdis,\n", + " ndmin,\n", + " ndmax,\n", + " radius,\n", + " ktype,\n", + " skmean,\n", + " vario,\n", + "):\n", + " UNEST = -999.\n", + " EPSLON = 1.0e-10\n", + " # VERSION = 2.907\n", + " # first = True\n", + " PMX = 9999.0 \n", + " MAXSAM = ndmax + 1\n", + " MAXDIS = nxdis * nydis * nzdis\n", + " MAXKD = MAXSAM + 1\n", + " MAXKRG = MAXKD * MAXKD\n", + "\n", + "# Allocate the needed memory: \n", + " xdb = np.zeros(MAXDIS)\n", + " ydb = np.zeros(MAXDIS)\n", + " zdb = np.zeros(MAXDIS)\n", + " xa = np.zeros(MAXSAM)\n", + " ya = np.zeros(MAXSAM)\n", + " za = np.zeros(MAXSAM)\n", + " vra = np.zeros(MAXSAM)\n", + " dist = np.zeros(MAXSAM)\n", + " nums = np.zeros(MAXSAM)\n", + " r = np.zeros(MAXKD)\n", + " rr = np.zeros(MAXKD)\n", + " s = np.zeros(MAXKD)\n", + " a = np.zeros(MAXKRG)\n", + " kmap = np.zeros((nx,ny,nz))\n", + " vmap = np.zeros((nx,ny,nz))\n", + "\n", + " # trim values outside tmin and tmax\n", + " df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] \n", + " assert df[vcol].max() <= tmax and df[vcol].min() >= tmin\n", + " nd = len(df_extract)\n", + " ndmax = min(ndmax,nd)\n", + " x = df_extract[xcol].values\n", + " y = df_extract[ycol].values\n", + " z = df_extract[zcol].values\n", + " vr = df_extract[vcol].values\n", + "\n", + "\n", + " # Allocate the needed memory: \n", + " xdb = np.zeros(MAXDIS)\n", + " ydb = np.zeros(MAXDIS)\n", + " zdb = np.zeros(MAXDIS)\n", + " xa = np.zeros(MAXSAM)\n", + " ya = np.zeros(MAXSAM)\n", + " za = np.zeros(MAXSAM)\n", + " vra = np.zeros(MAXSAM)\n", + " dist = np.zeros(MAXSAM)\n", + " nums = np.zeros(MAXSAM)\n", + " r = np.zeros(MAXKD)\n", + " rr = np.zeros(MAXKD)\n", + " s = np.zeros(MAXKD)\n", + " a = np.zeros(MAXKRG)\n", + " kmap = np.zeros((nx,ny,nz))\n", + " vmap = np.zeros((nx,ny,nz)) \n", + "\n", + "# set up tree for nearest neighbor search\n", + " dp = list((z[i], y[i], x[i]) for i in range(0,nd))\n", + " data_locs = np.column_stack((z, y, x))\n", + " tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error?\n", + "\n", + "# Summary statistics for the data after trimming\n", + " avg = vr.mean()\n", + " stdev = vr.std()\n", + " ss = stdev**2.0 # variance\n", + " vrmin = vr.min()\n", + " vrmax = vr.max()\n", + "\n", + "# load the variogram\n", + " nst = vario['nst'] # num structures\n", + " cc = np.zeros(nst) # variance contribution\n", + " aa = np.zeros(nst) # major range\n", + " it = np.zeros(nst) # type vario structure\n", + " aa = np.zeros(nst) # major range\n", + " ang_azi = np.zeros(nst) # azimuth\n", + " ang_dip = np.zeros(nst) # dip\n", + " anis = np.zeros(nst) # anistropy ratio w/ minor range\n", + " anis_v = np.zeros(nst) # anistropy ratio w/ vertical range\n", + " \n", + " # only works w/ nst == 1 or nst == 2\n", + " print(\"vario\")\n", + " print(vario)\n", + " c0 = vario['nug'] \n", + " it[0] = vario['it1'] \n", + " cc[0] = vario['cc1'] \n", + " ang_azi[0] = vario['azi1']\n", + " ang_dip[0] = vario['dip1']\n", + " aa[0] = vario['hmax']\n", + " anis[0] = vario['hmed1']/vario['hmax1']\n", + " anis_v[0] = vario['hmin1']/vario['hmax1']\n", + " if nst == 2:\n", + " cc[1] = vario['cc2']\n", + " it[1] = vario['it2']\n", + " ang_azi[1] = vario['azi2']\n", + " ang_dip[1] = vario['dip1']\n", + " aa[1] = vario['hmax']\n", + " anis[1] = vario['hmed2']/vario['hmax2']\n", + " anis_v[1] = vario['hmin2']/vario['hmax2']\n", + "\n", + " #rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX)\n", + " #cbb = maxcov\n", + "\n", + " # TODO Set up discretization points per block (bock kriging)\n", + "\n", + " # need to set up numPy cube grid\n", + " cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') \n", + " print(cubeGrid)\n", + " # main loop over all points:\n", + " nk = 0\n", + " ak = 0.0\n", + " vk = 0.0\n", + " for iz in range(0,nz):\n", + " zloc = zmn + (iz-0)*zsiz \n", + " for iy in range(0,ny):\n", + " yloc = ymn + (iy-0)*ysiz\n", + " for ix in range (0,nz): #triple nested for loop that loops through points and adds to matrices\n", + " xloc = xmn + (ix-0)*xsiz\n", + " current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell\n", + "\n", + "\n", + " # Find the nearest samples within each octant: First initialize\n", + "# the counter arrays:\n", + " na = -1 # accounting for 0 as first index\n", + " dist.fill(1.0e+20)\n", + " nums.fill(-1)\n", + " dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search\n", + " # remove any data outside search radius\n", + " na = len(dist)\n", + " nums = nums[dist UNEST:\n", + " nk = nk + 1\n", + " ak = ak + est\n", + " vk = vk + est*est\n", + "# TODO Set up kriging matrices:\n", + "\n", + " return kmap, vmap \n" + ] + }, + { + "cell_type": "code", + "execution_count": 100, + "metadata": {}, + "outputs": [], + "source": [ + "skmean_por = 0.10; skmean_perm = 65.0 # simple kriging mean (used if simple kriging is selected below)\n", + "ktype = 1 # kriging type, 0 - simple, 1 - ordinary\n", + "radius = 100 # search radius for neighbouring data\n", + "nxdis = 1; nydis = 1 ; nzdis = 1 # number of grid discretizations for block kriging (not tested)\n", + "ndmin = 0; ndmax = 10 # minimum and maximum data for an estimate\n", + "tmin = 0.0 # minimum property value\n", + "xsiz = 10; ysiz = 10 ; zsiz = 10 # cell size\n", + "nx = 100; ny = 100; nz = 100 # number of cells\n", + "xmn = 5; ymn = 5 ; zmn = 5 # grid origin, location center of lower left cell\n", + "por_vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0,hmaj1=100,hmin1=100) # porosity variogram\n", + "por_vario[\"dip1\"] = 0\n", + "por_vario[\"hmax\"] = 0\n", + "por_vario[\"hmed1\"] = 0\n", + "por_vario[\"hmax1\"] = 1\n" + ] + }, + { + "cell_type": "code", + "execution_count": 101, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "vario\n", + "{'nug': 0.0, 'nst': 1, 'it1': 1, 'cc1': 1.0, 'azi1': 0, 'hmaj1': 100, 'hmin1': 100, 'it2': 1, 'cc2': 0, 'azi2': 0, 'hmaj2': 0, 'hmin2': 0, 'dip1': 0, 'hmax': 0, 'hmed1': 0, 'hmax1': 1}\n", + "[[[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]\n", + "\n", + " ...\n", + "\n", + " [[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]\n", + "\n", + " [[0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " ...\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]\n", + " [0. 0. 0. ... 0. 0. 0.]]]\n" + ] + }, + { + "ename": "NameError", + "evalue": "name 'geostats' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mdf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"completions_production.csv\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity', tmin, tmax, nx,\n\u001b[1;32m----> 3\u001b[1;33m xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m\u001b[0m in \u001b[0;36mkt3d\u001b[1;34m(df, xcol, ycol, zcol, vcol, tmin, tmax, nx, xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean, vario)\u001b[0m\n\u001b[0;32m 176\u001b[0m \u001b[1;31m# Handle the situation of only one sample:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 177\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mna\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# accounting for min index of 0 - one sample case na = 0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 178\u001b[1;33m \u001b[0mcb1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgeostats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcova3\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mya\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mza\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mya\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mza\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnst\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mc0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mPMX\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0maa\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mit\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mang\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0manis\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mrotmat\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmaxcov\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 179\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 180\u001b[0m \u001b[0mxx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mxloc\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mNameError\u001b[0m: name 'geostats' is not defined" + ] + } + ], + "source": [ + "df = pd.read_csv(\"completions_production.csv\")\n", + "por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity', tmin, tmax, nx,\n", + " xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n", + "print(por_kmap)\n", + "print(por_vmap)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 24e2f5de6e5e615e1cce748afbcd389ada9cacf2 Mon Sep 17 00:00:00 2001 From: thingimajig13 <54867822+thingimajig13@users.noreply.github.com> Date: Tue, 5 May 2020 16:08:57 -0500 Subject: [PATCH 12/15] Delete testing.ipynb --- geostatspy/testing.ipynb | 471 --------------------------------------- 1 file changed, 471 deletions(-) delete mode 100644 geostatspy/testing.ipynb diff --git a/geostatspy/testing.ipynb b/geostatspy/testing.ipynb deleted file mode 100644 index fc3e27a..0000000 --- a/geostatspy/testing.ipynb +++ /dev/null @@ -1,471 +0,0 @@ -{ - "cells": [ - { - "cell_type": "code", - "execution_count": 88, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "testing the kt3d function using some data from Geostats Guy\n" - ] - } - ], - "source": [ - "print(\"testing the kt3d function using some data from Geostats Guy\")" - ] - }, - { - "cell_type": "code", - "execution_count": 98, - "metadata": {}, - "outputs": [ - { - "ename": "IndentationError", - "evalue": "unindent does not match any outer indentation level (geostats.py, line 364)", - "output_type": "error", - "traceback": [ - "Traceback \u001b[1;36m(most recent call last)\u001b[0m:\n", - " File \u001b[0;32m\"C:\\Users\\meena\\Anaconda3\\lib\\site-packages\\IPython\\core\\interactiveshell.py\"\u001b[0m, line \u001b[0;32m3325\u001b[0m, in \u001b[0;35mrun_code\u001b[0m\n exec(code_obj, self.user_global_ns, self.user_ns)\n", - "\u001b[1;36m File \u001b[1;32m\"\"\u001b[1;36m, line \u001b[1;32m3\u001b[1;36m, in \u001b[1;35m\u001b[1;36m\u001b[0m\n\u001b[1;33m import geostatspy.geostats as geostats\u001b[0m\n", - "\u001b[1;36m File \u001b[1;32m\"C:\\Users\\meena\\AppData\\Roaming\\Python\\Python37\\site-packages\\geostatspy\\geostats.py\"\u001b[1;36m, line \u001b[1;32m364\u001b[0m\n\u001b[1;33m EPSLON=1.0e-20\u001b[0m\n\u001b[1;37m ^\u001b[0m\n\u001b[1;31mIndentationError\u001b[0m\u001b[1;31m:\u001b[0m unindent does not match any outer indentation level\n" - ] - } - ], - "source": [ - "import geostatspy.GSLIB as GSLIB \n", - "# GSLIB utilies, visualization and wrapper\n", - "import geostatspy.geostats as geostats\n", - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "import scipy.spatial as sp # for fast nearest neighbor search" - ] - }, - { - "cell_type": "code", - "execution_count": 99, - "metadata": {}, - "outputs": [], - "source": [ - "\n", - "\n", - "\n", - "# *** all parameters have been copied over form kb2d with respective z parameter added if necessary\n", - "# we can probably delete a couple of these if we don't do block kriging\n", - "def kt3d (\n", - " df,\n", - " xcol,\n", - " ycol,\n", - " zcol,\n", - " vcol,\n", - " tmin,\n", - " tmax,\n", - " nx,\n", - " xmn,\n", - " xsiz,\n", - " ny,\n", - " ymn,\n", - " ysiz,\n", - " nz,\n", - " zmn,\n", - " zsiz,\n", - " nxdis,\n", - " nydis,\n", - " nzdis,\n", - " ndmin,\n", - " ndmax,\n", - " radius,\n", - " ktype,\n", - " skmean,\n", - " vario,\n", - "):\n", - " UNEST = -999.\n", - " EPSLON = 1.0e-10\n", - " # VERSION = 2.907\n", - " # first = True\n", - " PMX = 9999.0 \n", - " MAXSAM = ndmax + 1\n", - " MAXDIS = nxdis * nydis * nzdis\n", - " MAXKD = MAXSAM + 1\n", - " MAXKRG = MAXKD * MAXKD\n", - "\n", - "# Allocate the needed memory: \n", - " xdb = np.zeros(MAXDIS)\n", - " ydb = np.zeros(MAXDIS)\n", - " zdb = np.zeros(MAXDIS)\n", - " xa = np.zeros(MAXSAM)\n", - " ya = np.zeros(MAXSAM)\n", - " za = np.zeros(MAXSAM)\n", - " vra = np.zeros(MAXSAM)\n", - " dist = np.zeros(MAXSAM)\n", - " nums = np.zeros(MAXSAM)\n", - " r = np.zeros(MAXKD)\n", - " rr = np.zeros(MAXKD)\n", - " s = np.zeros(MAXKD)\n", - " a = np.zeros(MAXKRG)\n", - " kmap = np.zeros((nx,ny,nz))\n", - " vmap = np.zeros((nx,ny,nz))\n", - "\n", - " # trim values outside tmin and tmax\n", - " df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] \n", - " assert df[vcol].max() <= tmax and df[vcol].min() >= tmin\n", - " nd = len(df_extract)\n", - " ndmax = min(ndmax,nd)\n", - " x = df_extract[xcol].values\n", - " y = df_extract[ycol].values\n", - " z = df_extract[zcol].values\n", - " vr = df_extract[vcol].values\n", - "\n", - "\n", - " # Allocate the needed memory: \n", - " xdb = np.zeros(MAXDIS)\n", - " ydb = np.zeros(MAXDIS)\n", - " zdb = np.zeros(MAXDIS)\n", - " xa = np.zeros(MAXSAM)\n", - " ya = np.zeros(MAXSAM)\n", - " za = np.zeros(MAXSAM)\n", - " vra = np.zeros(MAXSAM)\n", - " dist = np.zeros(MAXSAM)\n", - " nums = np.zeros(MAXSAM)\n", - " r = np.zeros(MAXKD)\n", - " rr = np.zeros(MAXKD)\n", - " s = np.zeros(MAXKD)\n", - " a = np.zeros(MAXKRG)\n", - " kmap = np.zeros((nx,ny,nz))\n", - " vmap = np.zeros((nx,ny,nz)) \n", - "\n", - "# set up tree for nearest neighbor search\n", - " dp = list((z[i], y[i], x[i]) for i in range(0,nd))\n", - " data_locs = np.column_stack((z, y, x))\n", - " tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error?\n", - "\n", - "# Summary statistics for the data after trimming\n", - " avg = vr.mean()\n", - " stdev = vr.std()\n", - " ss = stdev**2.0 # variance\n", - " vrmin = vr.min()\n", - " vrmax = vr.max()\n", - "\n", - "# load the variogram\n", - " nst = vario['nst'] # num structures\n", - " cc = np.zeros(nst) # variance contribution\n", - " aa = np.zeros(nst) # major range\n", - " it = np.zeros(nst) # type vario structure\n", - " aa = np.zeros(nst) # major range\n", - " ang_azi = np.zeros(nst) # azimuth\n", - " ang_dip = np.zeros(nst) # dip\n", - " anis = np.zeros(nst) # anistropy ratio w/ minor range\n", - " anis_v = np.zeros(nst) # anistropy ratio w/ vertical range\n", - " \n", - " # only works w/ nst == 1 or nst == 2\n", - " print(\"vario\")\n", - " print(vario)\n", - " c0 = vario['nug'] \n", - " it[0] = vario['it1'] \n", - " cc[0] = vario['cc1'] \n", - " ang_azi[0] = vario['azi1']\n", - " ang_dip[0] = vario['dip1']\n", - " aa[0] = vario['hmax']\n", - " anis[0] = vario['hmed1']/vario['hmax1']\n", - " anis_v[0] = vario['hmin1']/vario['hmax1']\n", - " if nst == 2:\n", - " cc[1] = vario['cc2']\n", - " it[1] = vario['it2']\n", - " ang_azi[1] = vario['azi2']\n", - " ang_dip[1] = vario['dip1']\n", - " aa[1] = vario['hmax']\n", - " anis[1] = vario['hmed2']/vario['hmax2']\n", - " anis_v[1] = vario['hmin2']/vario['hmax2']\n", - "\n", - " #rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX)\n", - " #cbb = maxcov\n", - "\n", - " # TODO Set up discretization points per block (bock kriging)\n", - "\n", - " # need to set up numPy cube grid\n", - " cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') \n", - " print(cubeGrid)\n", - " # main loop over all points:\n", - " nk = 0\n", - " ak = 0.0\n", - " vk = 0.0\n", - " for iz in range(0,nz):\n", - " zloc = zmn + (iz-0)*zsiz \n", - " for iy in range(0,ny):\n", - " yloc = ymn + (iy-0)*ysiz\n", - " for ix in range (0,nz): #triple nested for loop that loops through points and adds to matrices\n", - " xloc = xmn + (ix-0)*xsiz\n", - " current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell\n", - "\n", - "\n", - " # Find the nearest samples within each octant: First initialize\n", - "# the counter arrays:\n", - " na = -1 # accounting for 0 as first index\n", - " dist.fill(1.0e+20)\n", - " nums.fill(-1)\n", - " dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search\n", - " # remove any data outside search radius\n", - " na = len(dist)\n", - " nums = nums[dist UNEST:\n", - " nk = nk + 1\n", - " ak = ak + est\n", - " vk = vk + est*est\n", - "# TODO Set up kriging matrices:\n", - "\n", - " return kmap, vmap \n" - ] - }, - { - "cell_type": "code", - "execution_count": 100, - "metadata": {}, - "outputs": [], - "source": [ - "skmean_por = 0.10; skmean_perm = 65.0 # simple kriging mean (used if simple kriging is selected below)\n", - "ktype = 1 # kriging type, 0 - simple, 1 - ordinary\n", - "radius = 100 # search radius for neighbouring data\n", - "nxdis = 1; nydis = 1 ; nzdis = 1 # number of grid discretizations for block kriging (not tested)\n", - "ndmin = 0; ndmax = 10 # minimum and maximum data for an estimate\n", - "tmin = 0.0 # minimum property value\n", - "xsiz = 10; ysiz = 10 ; zsiz = 10 # cell size\n", - "nx = 100; ny = 100; nz = 100 # number of cells\n", - "xmn = 5; ymn = 5 ; zmn = 5 # grid origin, location center of lower left cell\n", - "por_vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0,hmaj1=100,hmin1=100) # porosity variogram\n", - "por_vario[\"dip1\"] = 0\n", - "por_vario[\"hmax\"] = 0\n", - "por_vario[\"hmed1\"] = 0\n", - "por_vario[\"hmax1\"] = 1\n" - ] - }, - { - "cell_type": "code", - "execution_count": 101, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "vario\n", - "{'nug': 0.0, 'nst': 1, 'it1': 1, 'cc1': 1.0, 'azi1': 0, 'hmaj1': 100, 'hmin1': 100, 'it2': 1, 'cc2': 0, 'azi2': 0, 'hmaj2': 0, 'hmin2': 0, 'dip1': 0, 'hmax': 0, 'hmed1': 0, 'hmax1': 1}\n", - "[[[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]\n", - "\n", - " [[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]\n", - "\n", - " [[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]\n", - "\n", - " ...\n", - "\n", - " [[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]\n", - "\n", - " [[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]\n", - "\n", - " [[0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " ...\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]\n", - " [0. 0. 0. ... 0. 0. 0.]]]\n" - ] - }, - { - "ename": "NameError", - "evalue": "name 'geostats' is not defined", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[0mdf\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mread_csv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"completions_production.csv\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 2\u001b[0m por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity', tmin, tmax, nx,\n\u001b[1;32m----> 3\u001b[1;33m xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n\u001b[0m\u001b[0;32m 4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m\u001b[0m in \u001b[0;36mkt3d\u001b[1;34m(df, xcol, ycol, zcol, vcol, tmin, tmax, nx, xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean, vario)\u001b[0m\n\u001b[0;32m 176\u001b[0m \u001b[1;31m# Handle the situation of only one sample:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 177\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mna\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m:\u001b[0m \u001b[1;31m# accounting for min index of 0 - one sample case na = 0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 178\u001b[1;33m \u001b[0mcb1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgeostats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcova3\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mya\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mza\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mya\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mza\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnst\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mc0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mPMX\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0maa\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mit\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mang\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0manis\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mrotmat\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mmaxcov\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 179\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 180\u001b[0m \u001b[0mxx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mxa\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mxloc\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mNameError\u001b[0m: name 'geostats' is not defined" - ] - } - ], - "source": [ - "df = pd.read_csv(\"completions_production.csv\")\n", - "por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity', tmin, tmax, nx,\n", - " xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n", - "print(por_kmap)\n", - "print(por_vmap)\n", - "\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.7.3" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From 25bfbcb276f73a09decdb1aabd504debceb5d911 Mon Sep 17 00:00:00 2001 From: thingimajig13 <54867822+thingimajig13@users.noreply.github.com> Date: Tue, 5 May 2020 16:09:26 -0500 Subject: [PATCH 13/15] Add files via upload --- geostatspy/testing.ipynb | 389 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 389 insertions(+) create mode 100644 geostatspy/testing.ipynb diff --git a/geostatspy/testing.ipynb b/geostatspy/testing.ipynb new file mode 100644 index 0000000..9844bc2 --- /dev/null +++ b/geostatspy/testing.ipynb @@ -0,0 +1,389 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 102, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "testing the kt3d function using some data from Geostats Guy\n" + ] + } + ], + "source": [ + "print(\"testing the kt3d function using some data from Geostats Guy\")" + ] + }, + { + "cell_type": "code", + "execution_count": 103, + "metadata": {}, + "outputs": [], + "source": [ + "import geostatspy.GSLIB as GSLIB \n", + "# GSLIB utilies, visualization and wrapper\n", + "import geostatspy.geostats as geostats\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import scipy.spatial as sp # for fast nearest neighbor search" + ] + }, + { + "cell_type": "code", + "execution_count": 104, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "\n", + "\n", + "# *** all parameters have been copied over form kb2d with respective z parameter added if necessary\n", + "# we can probably delete a couple of these if we don't do block kriging\n", + "def kt3d (\n", + " df,\n", + " xcol,\n", + " ycol,\n", + " zcol,\n", + " vcol,\n", + " tmin,\n", + " tmax,\n", + " nx,\n", + " xmn,\n", + " xsiz,\n", + " ny,\n", + " ymn,\n", + " ysiz,\n", + " nz,\n", + " zmn,\n", + " zsiz,\n", + " nxdis,\n", + " nydis,\n", + " nzdis,\n", + " ndmin,\n", + " ndmax,\n", + " radius,\n", + " ktype,\n", + " skmean,\n", + " vario,\n", + "):\n", + " UNEST = -999.\n", + " EPSLON = 1.0e-10\n", + " # VERSION = 2.907\n", + " # first = True\n", + " PMX = 9999.0 \n", + " MAXSAM = ndmax + 1\n", + " MAXDIS = nxdis * nydis * nzdis\n", + " MAXKD = MAXSAM + 1\n", + " MAXKRG = MAXKD * MAXKD\n", + "\n", + "# Allocate the needed memory: \n", + " xdb = np.zeros(MAXDIS)\n", + " ydb = np.zeros(MAXDIS)\n", + " zdb = np.zeros(MAXDIS)\n", + " xa = np.zeros(MAXSAM)\n", + " ya = np.zeros(MAXSAM)\n", + " za = np.zeros(MAXSAM)\n", + " vra = np.zeros(MAXSAM)\n", + " dist = np.zeros(MAXSAM)\n", + " nums = np.zeros(MAXSAM)\n", + " r = np.zeros(MAXKD)\n", + " rr = np.zeros(MAXKD)\n", + " s = np.zeros(MAXKD)\n", + " a = np.zeros(MAXKRG)\n", + " kmap = np.zeros((nx,ny,nz))\n", + " vmap = np.zeros((nx,ny,nz))\n", + "\n", + " # trim values outside tmin and tmax\n", + " df_extract = df.loc[(df[vcol] >= tmin) & (df[vcol] <= tmax)] \n", + " assert df[vcol].max() <= tmax and df[vcol].min() >= tmin\n", + " nd = len(df_extract)\n", + " ndmax = min(ndmax,nd)\n", + " x = df_extract[xcol].values\n", + " y = df_extract[ycol].values\n", + " z = df_extract[zcol].values\n", + " vr = df_extract[vcol].values\n", + "\n", + "\n", + " # Allocate the needed memory: \n", + " xdb = np.zeros(MAXDIS)\n", + " ydb = np.zeros(MAXDIS)\n", + " zdb = np.zeros(MAXDIS)\n", + " xa = np.zeros(MAXSAM)\n", + " ya = np.zeros(MAXSAM)\n", + " za = np.zeros(MAXSAM)\n", + " vra = np.zeros(MAXSAM)\n", + " dist = np.zeros(MAXSAM)\n", + " nums = np.zeros(MAXSAM)\n", + " r = np.zeros(MAXKD)\n", + " rr = np.zeros(MAXKD)\n", + " s = np.zeros(MAXKD)\n", + " a = np.zeros(MAXKRG)\n", + " kmap = np.zeros((nx,ny,nz))\n", + " vmap = np.zeros((nx,ny,nz)) \n", + "\n", + "# set up tree for nearest neighbor search\n", + " dp = list((z[i], y[i], x[i]) for i in range(0,nd))\n", + " data_locs = np.column_stack((z, y, x))\n", + " tree = sp.cKDTree(data_locs, leafsize=16, compact_nodes=True, copy_data=False, balanced_tree=True) #why this error?\n", + "\n", + "# Summary statistics for the data after trimming\n", + " avg = vr.mean()\n", + " stdev = vr.std()\n", + " ss = stdev**2.0 # variance\n", + " vrmin = vr.min()\n", + " vrmax = vr.max()\n", + "\n", + "# load the variogram\n", + " nst = vario['nst'] # num structures\n", + " cc = np.zeros(nst) # variance contribution\n", + " aa = np.zeros(nst) # major range\n", + " it = np.zeros(nst) # type vario structure\n", + " aa = np.zeros(nst) # major range\n", + " ang_azi = np.zeros(nst) # azimuth\n", + " ang_dip = np.zeros(nst) # dip\n", + " anis = np.zeros(nst) # anistropy ratio w/ minor range\n", + " anis_v = np.zeros(nst) # anistropy ratio w/ vertical range\n", + " \n", + " # only works w/ nst == 1 or nst == 2\n", + " print(\"vario\")\n", + " print(vario)\n", + " c0 = vario['nug'] \n", + " it[0] = vario['it1'] \n", + " cc[0] = vario['cc1'] \n", + " ang_azi[0] = vario['azi1']\n", + " ang_dip[0] = vario['dip1']\n", + " aa[0] = vario['hmax']\n", + " anis[0] = vario['hmed1']/vario['hmax1']\n", + " anis_v[0] = vario['hmin1']/vario['hmax1']\n", + " if nst == 2:\n", + " cc[1] = vario['cc2']\n", + " it[1] = vario['it2']\n", + " ang_azi[1] = vario['azi2']\n", + " ang_dip[1] = vario['dip1']\n", + " aa[1] = vario['hmax']\n", + " anis[1] = vario['hmed2']/vario['hmax2']\n", + " anis_v[1] = vario['hmin2']/vario['hmax2']\n", + "\n", + " #rotmat, maxcov = setup_rotmat_3D(c0, nst, it, cc, ang_azi, ang_dip, PMX)\n", + " #cbb = maxcov\n", + "\n", + " # TODO Set up discretization points per block (bock kriging)\n", + "\n", + " # need to set up numPy cube grid\n", + " cubeGrid = np.ndarray(shape=(nx,ny,nz), dtype=float, order='C') \n", + " print(cubeGrid)\n", + " # main loop over all points:\n", + " nk = 0\n", + " ak = 0.0\n", + " vk = 0.0\n", + " for iz in range(0,nz):\n", + " zloc = zmn + (iz-0)*zsiz \n", + " for iy in range(0,ny):\n", + " yloc = ymn + (iy-0)*ysiz\n", + " for ix in range (0,nz): #triple nested for loop that loops through points and adds to matrices\n", + " xloc = xmn + (ix-0)*xsiz\n", + " current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell\n", + "\n", + "\n", + " # Find the nearest samples within each octant: First initialize\n", + "# the counter arrays:\n", + " na = -1 # accounting for 0 as first index\n", + " dist.fill(1.0e+20)\n", + " nums.fill(-1)\n", + " dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search\n", + " # remove any data outside search radius\n", + " na = len(dist)\n", + " nums = nums[dist UNEST:\n", + " nk = nk + 1\n", + " ak = ak + est\n", + " vk = vk + est*est\n", + "# TODO Set up kriging matrices:\n", + "\n", + " return kmap, vmap \n" + ] + }, + { + "cell_type": "code", + "execution_count": 105, + "metadata": {}, + "outputs": [], + "source": [ + "skmean_por = 0.10; skmean_perm = 65.0 # simple kriging mean (used if simple kriging is selected below)\n", + "ktype = 1 # kriging type, 0 - simple, 1 - ordinary\n", + "radius = 100 # search radius for neighbouring data\n", + "nxdis = 1; nydis = 1 ; nzdis = 1 # number of grid discretizations for block kriging (not tested)\n", + "ndmin = 0; ndmax = 10 # minimum and maximum data for an estimate\n", + "tmin = 0.0 # minimum property value\n", + "xsiz = 10; ysiz = 10 ; zsiz = 10 # cell size\n", + "nx = 100; ny = 100; nz = 100 # number of cells\n", + "xmn = 5; ymn = 5 ; zmn = 5 # grid origin, location center of lower left cell\n", + "por_vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0,hmaj1=100,hmin1=100) # porosity variogram\n", + "por_vario[\"dip1\"] = 0\n", + "por_vario[\"hmax\"] = 0\n", + "por_vario[\"hmed1\"] = 0\n", + "por_vario[\"hmax1\"] = 1\n" + ] + }, + { + "cell_type": "code", + "execution_count": 106, + "metadata": {}, + "outputs": [], + "source": [ + "df = pd.read_csv(\"completions_production.csv\")\n", + "por_kmap, por_vmap = kt3d ( df, 'X', 'Y', 'Z', 'Porosity', tmin, tmax, nx,\n", + " xmn, xsiz, ny, ymn, ysiz, nz, zmn, zsiz, nxdis, nydis, nzdis, ndmin, ndmax, radius, ktype, skmean_por, por_vario)\n", + "print(por_kmap)\n", + "print(por_vmap)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 57233b83745decb701788849a88bed701e90bc52 Mon Sep 17 00:00:00 2001 From: fang-helen <53734084+fang-helen@users.noreply.github.com> Date: Tue, 5 May 2020 16:37:23 -0500 Subject: [PATCH 14/15] small fixes --- geostatspy/geostats.py | 25 +++++++++++++------------ 1 file changed, 13 insertions(+), 12 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index 9c680ca..a687aba 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2245,10 +2245,10 @@ def kt3d ( vra = np.zeros(MAXSAM) dist = np.zeros(MAXSAM) nums = np.zeros(MAXSAM) - r = np.zeros(MAXKD) - rr = np.zeros(MAXKD) - s = np.zeros(MAXKD) - a = np.zeros(MAXKRG) + # r = np.zeros(MAXKD) + # rr = np.zeros(MAXKD) + # s = np.zeros(MAXKD) + # a = np.zeros(MAXKRG) kmap = np.zeros((nx,ny,nz)) vmap = np.zeros((nx,ny,nz)) @@ -2269,12 +2269,12 @@ def kt3d ( ya = np.zeros(MAXSAM) za = np.zeros(MAXSAM) vra = np.zeros(MAXSAM) - dist = np.zeros(MAXSAM) - nums = np.zeros(MAXSAM) - r = np.zeros(MAXKD) - rr = np.zeros(MAXKD) - s = np.zeros(MAXKD) - a = np.zeros(MAXKRG) + # dist = np.zeros(MAXSAM) # these allocations shouldn't be necessary + # nums = np.zeros(MAXSAM) + # r = np.zeros(MAXKD) + # rr = np.zeros(MAXKD) + # s = np.zeros(MAXKD) + # a = np.zeros(MAXKRG) kmap = np.zeros((nx,ny,nz)) vmap = np.zeros((nx,ny,nz)) @@ -2343,8 +2343,8 @@ def kt3d ( # Find the nearest samples within each octant: First initializ the counter arrays: na = -1 # accounting for 0 as first index - dist.fill(1.0e+20) - nums.fill(-1) + # dist.fill(1.0e+20) + # nums.fill(-1) dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search # remove any data outside search radius na = len(dist) @@ -2403,6 +2403,7 @@ def kt3d ( # Set up kriging matrices: a = np.zeroes((na,na)) r = np.zeroes(na) + rr = np.zeroes(na) # Establish Left Hand Side Covariance Matrix: From 6bdf714b6b3952df9b3a74fd665a0aeb0fb009dd Mon Sep 17 00:00:00 2001 From: fang-helen Date: Sun, 10 May 2020 14:09:23 -0500 Subject: [PATCH 15/15] cleaning up --- geostatspy/geostats.py | 67 ++++++++++++------------------------------ 1 file changed, 18 insertions(+), 49 deletions(-) diff --git a/geostatspy/geostats.py b/geostatspy/geostats.py index a687aba..033feb1 100644 --- a/geostatspy/geostats.py +++ b/geostatspy/geostats.py @@ -2227,28 +2227,24 @@ def kt3d ( ): UNEST = -999. EPSLON = 1.0e-10 - # VERSION = 2.907 - # first = True PMX = 9999.0 MAXSAM = ndmax + 1 MAXDIS = nxdis * nydis * nzdis MAXKD = MAXSAM + 1 MAXKRG = MAXKD * MAXKD -# Allocate the needed memory: +# Allocate the needed memory: + + # used for block kriging displacements xdb = np.zeros(MAXDIS) ydb = np.zeros(MAXDIS) zdb = np.zeros(MAXDIS) + # temp arrays to hold nearest neighbor search results xa = np.zeros(MAXSAM) ya = np.zeros(MAXSAM) za = np.zeros(MAXSAM) vra = np.zeros(MAXSAM) - dist = np.zeros(MAXSAM) - nums = np.zeros(MAXSAM) - # r = np.zeros(MAXKD) - # rr = np.zeros(MAXKD) - # s = np.zeros(MAXKD) - # a = np.zeros(MAXKRG) + # result grid meshes for estimated values and estimation variance kmap = np.zeros((nx,ny,nz)) vmap = np.zeros((nx,ny,nz)) @@ -2261,23 +2257,6 @@ def kt3d ( z = df_extract[zcol].values vr = df_extract[vcol].values - # Allocate the needed memory: - xdb = np.zeros(MAXDIS) - ydb = np.zeros(MAXDIS) - zdb = np.zeros(MAXDIS) - xa = np.zeros(MAXSAM) - ya = np.zeros(MAXSAM) - za = np.zeros(MAXSAM) - vra = np.zeros(MAXSAM) - # dist = np.zeros(MAXSAM) # these allocations shouldn't be necessary - # nums = np.zeros(MAXSAM) - # r = np.zeros(MAXKD) - # rr = np.zeros(MAXKD) - # s = np.zeros(MAXKD) - # a = np.zeros(MAXKRG) - kmap = np.zeros((nx,ny,nz)) - vmap = np.zeros((nx,ny,nz)) - # set up tree for nearest neighbor search dp = list((z[i], y[i], x[i]) for i in range(0,nd)) data_locs = np.column_stack((z, y, x)) @@ -2340,11 +2319,8 @@ def kt3d ( xloc = xmn + (ix-0)*xsiz current_node = (zloc, yloc,xloc) # xloc, yloc, zloc centroid of current cell - - # Find the nearest samples within each octant: First initializ the counter arrays: + # Find the nearest samples within each octant: na = -1 # accounting for 0 as first index - # dist.fill(1.0e+20) - # nums.fill(-1) dist, nums = tree.query(current_node,ndmax) # use kd tree for fast nearest data search # remove any data outside search radius na = len(dist) @@ -2356,6 +2332,7 @@ def kt3d ( if na + 1 < ndmin: # accounting for min index of 0 est = UNEST estv = UNEST + # not enough samples, mark as "unestimated" print('UNEST at ' + str(ix) + ',' + str(iy) + ','+str(iz)) else: # Put coordinates and values of neighborhood samples into xa,ya,za,vra: @@ -2364,7 +2341,7 @@ def kt3d ( xa[ia] = x[jj] #why this error? ya[ia] = y[jj] za[ia] = z[jj] - vra[ia] = vr[jj] #what is vra?? + vra[ia] = vr[jj] # Handle the situation of only one sample: if na == 0: # accounting for min index of 0 - one sample case na = 0 @@ -2375,25 +2352,15 @@ def kt3d ( zz = za[0] - zloc # Establish Right Hand Side Covariance: - if ndb <= 1: - cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) - else: - cb = 0.0 - for i in range(0,ndb): - cb = cb + cova3(xx,yy,zz,xdb[i],ydb[i],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) - dx = xx - xdb(i) - dy = yy - ydb(i) - dz = zz - zdb(i) - if (dx*dx+dy*dy+dz*dz) < EPSLON: - cb = cb - c0 - cb = cb / real(ndb) + cb = cova3(xx,yy,zz,xdb[0],ydb[0],zdb[i],nst,c0,PMX,cc,aa,it,ang,anis,rotmat,maxcov) + if ktype == 0: s[0] = cb/cbb est = s[0]*vra[0] + (1.0-s[0])*skmean estv = cbb - s[0] * cb - else: - est = vra[0] - estv = cbb - 2.0*cb + cb1 + # else: + # est = vra[0] + # estv = cbb - 2.0*cb + cb1 else: # Solve the Kriging System with more than one sample: neq = na + ktype # accounting for first index of 0 @@ -2408,10 +2375,13 @@ def kt3d ( # Establish Left Hand Side Covariance Matrix: for j in range(0,na): - for i in range(0,na): #loops through and adds covariances into matrix + for i in range(0,na): + # add covariance into left-hand square matrix a[j][i] = cova3(xa[i],ya[i],za[i],xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) + # add covariance into right-hand column matrix r[j] = cova3(xloc,yloc,zloc,xa[j],ya[j],za[j],nst,c0,PMX,cc,aa,it,ang_azi,anis,rotmat,maxcov) + # solve the kriging system s = ksol_numpy(neq,a,r) ising = 0 @@ -2421,7 +2391,6 @@ def kt3d ( est = UNEST estv = UNEST else: - # Compute the estimate and the kriging variance: est = 0.0 estv = cbb @@ -2434,7 +2403,7 @@ def kt3d ( estv = estv - s[i]*rr[i] if ktype == 0: est = est + (1.0-sumw)*skmean - # add estimation variance and estimated value to map + # add estimation variance and estimated value to maps kmap[nz-iz-1, ny-iy-1,ix] = est vmap[nz-iz-1, ny-iy-1,ix] = estv if est > UNEST: