地理可视化绘制王炸工具！！推荐收藏~~

import numpy as np
import matplotlib.pyplot as plt

import pygimli as pg
from pygimli.viewer.mpl import drawMesh, drawModel
#Create coarse and fine mesh with data
def create_mesh_and_data(n):
    nc = np.linspace(-2.0, 0.0, n)
    mesh = pg.meshtools.createMesh2D(nc, nc)
    mcx = pg.x(mesh.cellCenter())
    mcy = pg.y(mesh.cellCenter())
    data = np.cos(1.5 * mcx) * np.sin(1.5 * mcy)
    return mesh, data

coarse, coarse_data = create_mesh_and_data(5)
fine, fine_data = create_mesh_and_data(20)
#Interpolate data to different meshes
def nearest_neighbor_interpolation(inmesh, indata, outmesh, nan=99.9):
    """ Nearest neighbor interpolation. """
    outdata = []
    for pos in outmesh.cellCenters():
        cell = inmesh.findCell(pos)
        if cell:
            outdata.append(indata[cell.id()])
        else:
            outdata.append(nan)
    return outdata


def linear_interpolation(inmesh, indata, outmesh):
    """ Linear interpolation using `pg.interpolate()` """
    outdata = pg.Vector()  # empty
    pg.interpolate(srcMesh=inmesh, inVec=indata,
                   destPos=outmesh.cellCenters(), outVec=outdata)

    # alternatively you can use the interpolation matrix
    outdata = inmesh.interpolationMatrix(outmesh.cellCenters()) * \
              pg.core.cellDataToPointData(inmesh, indata)
    return outdata
#Visualization
meshes = [coarse, fine]
datasets = [coarse_data, fine_data]
ints = [nearest_neighbor_interpolation,
        linear_interpolation]

fig, ax = plt.subplots(2, 2, figsize=(5, 5))

# Coarse data to fine mesh
drawModel(ax[0, 0], fine, ints[0](coarse, coarse_data, fine), showCbar=False)
drawMesh(ax[0, 0], fine)
drawModel(ax[0, 1], fine, ints[1](coarse, coarse_data, fine), showCbar=False)
drawMesh(ax[0, 1], fine)

# Fine data to coarse mesh
drawModel(ax[1, 0], coarse, ints[0](fine, fine_data, coarse), showCbar=False)
drawMesh(ax[1, 0], coarse)
drawModel(ax[1, 1], coarse, ints[1](fine, fine_data, coarse), showCbar=False)
drawMesh(ax[1, 1], coarse)

titles = ["Coarse to fine\nwith nearest neighbors",
          "Coarse to fine\nwith linear interpolation",
          "Fine to coarse\nwith nearest neighbors",
          "Fine to coarse\nwith linear interpolation"]

for a, title in zip(ax.flat, titles):
    a.set_title(title + "\n")

fig.tight_layout()
plt.show()

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
z = np.geomspace(1, 5, 5)-1
cylinder = pg.meshtools.extrudeMesh(circle, a=z)
pg.show(cylinder, cylinder.cellMarkers(), label="Cell markers")

import numpy as np
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import traveltime
from pygimli.viewer.pv import drawSensors

pyvista = pg.optImport("pyvista")
depth = 15
width = 30
plc = mt.createCube(size=[width, width, depth], pos=[0, 0, -depth/2], area=5)

n_sensors = 8
sensors = np.zeros((n_sensors, 3))
sensors[0, 0] = 15
sensors[0, 1] = -10
sensors[1:, 0] = -15
sensors[1:, 1] = np.linspace(-15, 15, n_sensors - 1)

for pos in sensors:
    plc.createNode(pos)
mesh = mt.createMesh(plc)
mesh.createSecondaryNodes(1)
vel = 300 + -pg.z(mesh.cellCenters()) * 100

if pyvista:
    label = pg.utils.unit("vel")
    pg.show(mesh, vel, label=label)

import numpy as np

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
# from pygimli.physics.ert import simulate as simulateERT
from pygimli.physics.ert import VESModelling, VESCModelling
# from pygimli.physics.ert import createERTData

scheme = ert.createData(pg.utils.grange(start=1, end=24, dx=1, n=10, log=True),
                        sounding=True)
plc = mt.createWorld(start=[-200, -100], end=[200, 0],
                     layers=[-10], area=[5.0, 500])
for s in scheme.sensors():
    plc.createNode(s + [0.0, -0.2])

# Now we can create our forward modeling mesh.
mesh = mt.createMesh(plc, quality=33)

pg.show(mesh, data=mesh.cellMarkers(), label='Marker', showMesh=True)

import numpy as np

import pygimli as pg
import pygimli.meshtools as mt
from pygimli.viewer.pv import drawStreamLines, drawSlice

plc = mt.createCube(size=[40, 20, 15], marker=1, boundaryMarker=0)
cube = mt.createCube(size=[15, 15, 8], marker=2, boundaryMarker=0)
geom = plc + cube

mesh = mt.createMesh(geom, area=4)

for bound in mesh.boundaries():
    x = bound.center().x()
    if x == mesh.xmin():
        bound.setMarker(1)
    elif x == mesh.xmax():
        bound.setMarker(2)

kMap = {1: 1e-4, 2: 1e-6}
kArray = pg.solver.parseMapToCellArray(list(kMap), mesh) # dict does not work
kArray = np.column_stack([kArray] * 3)

bc = {"Dirichlet": {1: 20.0, 2: 10.0}}


h = pg.solver.solveFiniteElements(mesh, kMap, bc=bc)
vel = -pg.solver.grad(mesh, h) * kArray

pg.show(mesh, h, label="Hydraulic head (m)")

ax, _ = pg.show(mesh, hold=True, alpha=0.3)
drawStreamLines(ax, mesh, vel, radius=.1, source_radius=10)
drawSlice(ax, mesh, normal=[0,1,0], data=pg.abs(vel), label="Absolute velocity")
ax.show()