Plot 3D Meshing¶

In [1]:
%%javascript
IPython.OutputArea.prototype._should_scroll = function(lines) {
return false;
}
In [2]:
import numpy as np     
import plotly.graph_objects as go
import gstlearn.plot3D as gop
import IPython
import os
import urllib.request
from numpy import pi, cos, sin
import gstlearn as gl

On the Sphere¶

Definition of the Meshing

In [3]:
gl.defineDefaultSpace(gl.ESpaceType.SN)
mesh = gl.MeshSphericalExt()
err = mesh.resetFromDb(None,None,triswitch = "-r4",verbose=False)

Display a white skin around the meshing

In [4]:
blank = gop.SurfaceOnMesh(mesh, opacity=1)
In [5]:
fig = go.Figure(data = [blank])
fig.update_scenes(xaxis_visible=False, yaxis_visible=False, zaxis_visible=False )
f = fig.show()

We overlay the meshing

In [6]:
meshing = gop.Meshing(mesh)
In [7]:
fig = go.Figure(data = [blank, meshing])
fig.update_scenes(xaxis_visible=False, yaxis_visible=False, zaxis_visible=False )
f = fig.show()

Drawing a polygon (we use the one containing the land boundaries)

In [8]:
url = 'https://soft.minesparis.psl.eu/gstlearn/data/boundaries/world.poly'
name, head = urllib.request.urlretrieve(url)
poly = gl.Polygons.createFromNF(name)
poly.display()

Polygons
--------
Number of Polygon Sets = 451
In [9]:
boundaries = gop.PolygonOnSphere(poly)
equator    = gop.Equator(width=5)
meridians  = gop.Meridians(angle=20,color='blue')
parallels  = gop.Parallels(angle=30,color='red')
pole       = gop.Pole()
poleaxis   = gop.PolarAxis()
In [10]:
fig = go.Figure(data = [blank,boundaries,equator,meridians,parallels,pole,poleaxis])
fig.update_scenes(xaxis_visible=False, yaxis_visible=False, zaxis_visible=False )
f = fig.show()

Representing a function on the skin of the earth together with other decoration. The function is the result of a non-conditional simulation performed with SPDE on the Sphere.

In [11]:
model = gl.Model.createFromParam(gl.ECov.BESSEL_K,range=1500,param=1)
S = gl.ShiftOpCs(mesh,model)
Q = gl.PrecisionOpCs(S,model.getCova(0))
result = Q.simulateOne()
In [12]:
simu = gop.SurfaceOnMesh(mesh, intensity=result, opacity=1)
In [13]:
fig = go.Figure(data = [simu,boundaries,equator,meridians,parallels,pole,poleaxis])
fig.update_scenes(xaxis_visible=False, yaxis_visible=False, zaxis_visible=False )
f = fig.show()

3D in General¶

We define the space dimension

In [14]:
ndim = 3
gl.defineDefaultSpace(gl.ESpaceType.RN, ndim)

Defining the output grid

In [15]:
nx = [61,81,61]
dx = [0.1, 0.1, 0.1]
x0 = [-3., -4., -6.]
grid = gl.DbGrid.create(nx=nx, dx=dx, x0=x0)
x = grid.getCoordinates(0)
y = grid.getCoordinates(1)
z = grid.getCoordinates(2)
val = x*x + y*y + z*z
grid.addColumns(val,"Data",gl.ELoc.Z)
glimits = grid.getRange("Data")

Defining a Data Set with Gradient and Tangent components

In [16]:
nech = 5
coormin = grid.getCoorMinimum()
coormax = grid.getCoorMaximum()
db = gl.Db.createFromBox(nech, coormin, coormax)
np.random.seed(123)

# Data
uid = db.addColumns(np.random.uniform(10.,20.,nech), "Data", gl.ELoc.Z)

# Gradient components
uid = db.addColumns(np.random.normal(0, 1, nech),"gx",gl.ELoc.G,0)
uid = db.addColumns(np.random.normal(0, 1, nech),"gy",gl.ELoc.G,1)
uid = db.addColumns(np.random.normal(0, 1, nech),"gz",gl.ELoc.G,2)

# Tangent components
uid = db.addColumns(np.random.normal(0, 1, nech),"tx",gl.ELoc.TGTE,0)
uid = db.addColumns(np.random.normal(0, 1, nech),"ty",gl.ELoc.TGTE,1)
uid = db.addColumns(np.random.normal(0, 1, nech),"tz",gl.ELoc.TGTE,2)

3-D Visualization

In [17]:
levels = [5., 25.]
surf1 = gop.IsoSurfaceOnDbGrid(grid, "Data", useSel=False, isomin=levels[0], isomax=levels[0])
surf2 = gop.IsoSurfaceOnDbGrid(grid, "Data", useSel=False, isomin=levels[1], isomax=levels[1])
point = gop.PointDb(db, size=5, name_color = "Data")
gradient = gop.GradientDb(db,size=0.5,colorscale='blues',sizemode='absolute')
tangent  = gop.TangentDb(db,size=0.5,colorscale='gray',sizemode='absolute')
In [18]:
fig = go.Figure(data = [ surf1, surf2, point, gradient, tangent])
f = fig.show()