前言

采用 Mines JTK 绘图的好处:

  • 符合 Geophysics 等主流杂志的图片要求
  • 自动设置图片的各种参数,并根据需求(论文,演讲)生成不同大小的 label 文字
  • 三维绘图能力突出,个人感觉比 MATLAB 要美观。

Mines JTK 基础配置

二维图形

import sys
import os

from java.awt import *
from java.io import *
from java.lang import *
from javax.swing import *
from java.util import *

# MINES_JTK_LIB 为 edu-mines-jtk-1.1.0.jar 文件的绝对位置,此例中将其设为系统环境变量。
sys.path += [os.environ.get('MINES_JTK_LIB')]

from edu.mines.jtk.awt import *
from edu.mines.jtk.dsp import *
from edu.mines.jtk.io import *
from edu.mines.jtk.interp import *
from edu.mines.jtk.mosaic import *
from edu.mines.jtk.util import *
from edu.mines.jtk.util.ArrayMath import *

pngDir = None
pngDir = "./png/"

seismicDir = "./data/"
fxfile = "vel2"
f1,f2 = 0,0
d1,d2 = 0.012192,0.012192
n1,n2 = 300,1290
s1 = Sampling(n1,d1,f1)
s2 = Sampling(n2,d2,f2)

def main(args):
  greyPlot()
  #colorPlot()

def greyPlot():
  fx = readImage(fxfile)
  # 绘制灰度图
  plot(fx,cmap=gray,cmin=1.500,cmax=4.500,cint=0.500,label="Velocity (km/s)",png="vel-grey")
  # 绘制彩色图 
  plot(fx,cmap=jet,cmin=1.500,cmax=4.500,cint=0.500,label="Velocity (km/s)",png="vel-jet")
#############################################################################
# graphics
gray = ColorMap.GRAY
jet = ColorMap.JET
def jetFill(alpha):
  return ColorMap.setAlpha(ColorMap.JET,alpha)
def jetFillExceptMin(alpha):
  a = fillfloat(alpha,256)
  a[0] = 0.0
  return ColorMap.setAlpha(ColorMap.JET,a)
def jetRamp(alpha):
  return ColorMap.setAlpha(ColorMap.HUE_BLUE_TO_RED,rampfloat(0.0,alpha/256,256))
def grayRamp(alpha):
  return ColorMap.setAlpha(ColorMap.GRAY,rampfloat(0.0,alpha/256,256))

def plot(f,g=None,v1=None,v2=None, cmap=None,cmin=None,cmax=None,cint=None,
        label=None,neareast=False,png=None):
  sp = SimplePlot(SimplePlot.Origin.UPPER_LEFT)
  sp.setVInterval(1.000)
  sp.setHInterval(2.000)
  sp.setHLabel("Distance (km)")
  sp.setVLabel("Depth (km)")
  pxv = sp.addPixels(s1,s2,f)
  sp.setVLimits(f1,f1+n1*d1)
  sp.setHLimits(f2,f2+n2*d2)
  if cmap:
    pxv.setColorModel(cmap)
  else:
    pxv.setColorModel(ColorMap.JET)
  #pxv.setInterpolation(PixelsView.Interpolation.NEAREST)
  if g:
    pxv.setClips(-1,1)
  else:
    if cmin and cmax:
      pxv.setClips(cmin,cmax)
  if g:
    pv = sp.addPixels(s1,s2,g)
    pv.setInterpolation(PixelsView.Interpolation.NEAREST)
    pv.setColorModel(cmap)
    if cmin and cmax:
      pv.setClips(cmin,cmax)
  cb = sp.addColorBar()
  if cint:
    cb.setInterval(cint)
  if label:
    cb.setLabel(label)
  if (v1 and v2):
    x1 = zerofloat(2)
    x2 = zerofloat(2)
    dx1,dx2 = 6, 14
    scale = 4
    for i2 in range(0,n2,dx2):
      for i1 in range(0,n1,dx1):
        x2[0] = (i2-v2[i2][i1]*scale)*d2+f2
        x2[1] = (i2+v2[i2][i1]*scale)*d2+f2
        x1[0] = (i1-v1[i2][i1]*scale)*d1+f1
        x1[1] = (i1+v1[i2][i1]*scale)*d1+f1
        pvu = sp.addPoints(x1,x2)
        pvu.setLineWidth(4)
        pvu.setLineColor(Color.CYAN)
  sp.plotPanel.setColorBarWidthMinimum(100)
  sp.setSize(1050,700) #for f3d
  sp.setFontSize(30)
  if pngDir and png:
    sp.paintToPng(300,3.333,pngDir+png+".png")


#############################################################################
# utilities
def gain(x):
  g = mul(x,x)
  ref = RecursiveExponentialFilter(100.0)
  ref.apply1(g,g)
  y = zerofloat(n1,n2)
  div(x,sqrt(g),y)
  return y

def readImage(name):
  fileName = seismicDir+name+".dat" # this data must be xdr_float format.
  n1,n2 = s1.count,s2.count
  image = zerofloat(n1,n2)
  ais = ArrayInputStream(fileName)
  ais.readFloats(image)
  ais.close()
  return image

#############################################################################
# Run the function main on the Swing thread
import sys
class _RunMain(Runnable):
  def __init__(self,main):
    self.main = main
  def run(self):
    self.main(sys.argv)
def run(main):
  SwingUtilities.invokeLater(_RunMain(main))
run(main)

vel-grey.png

vel-jet.png

参考资料

  1. https://github.com/xinwucwp/sos
  2. http://inside.mines.edu/~dhale/notebook.html#2009_07_05