p-v图是分析分子云中气体动力学的一种常用工具。简单讲就是沿某个特定方向分析谱线形态的变化，从而了解此方向上气体的运动规律。p-v图经常用来辅助寻找外流（高速）气体。
       今天写了个小程序可以画出一个正方形区域的过中心的横、竖、两对角共四个方向的p-v图。
       另附上python画图的颜色表以备不时之需（通过修改cmap=cm.spectral中"."后面的名称选择颜色表，http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps）。

===========================================================
#!/usr/bin/env python

import pyfits
from pylab import *
import sys

import math
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import mpl_toolkits
from matplotlib.patches import Ellipse

hdulist = pyfits.open('xxx.fits')
image = hdulist[0].data
nx = hdulist[0].header['naxis1']
ny = hdulist[0].header['naxis2']
nz = hdulist[0].header['naxis3']
crvalx = hdulist[0].header['crval1']
cdeltax = hdulist[0].header['cdelt1']
crpixx = hdulist[0].header['crpix1']
crvaly = hdulist[0].header['crval2']
cdeltay = hdulist[0].header['cdelt2']
crpixy = hdulist[0].header['crpix2']
crvalz = hdulist[0].header['crval3']
cdeltaz = hdulist[0].header['cdelt3']
crpixz = hdulist[0].header['crpix3']

x = np.arange(-crpixx*cdeltax+crvalx,(nx-1-crpixx)*cdeltax+crvalx,cdeltax)
y = np.arange(-crpixy*cdeltay+crvaly,(ny-1-crpixy)*cdeltay+crvaly,cdeltay)
z = np.arange(-crpixz*cdeltaz+crvalz,(nz-1-crpixz)*cdeltaz+crvalz,cdeltaz)
z = z/1000.0

# image[z,y,x]!!!! The first dimension is velocity channel!

print 'input'
print 'centerx,centery,xscale,yscale'
line = sys.stdin.readline()
[centerx,centery,xscale,yscale] = line.rstrip().split(',')
centerx=int(centerx)
centery=int(centery)
xscale=int(xscale)
yscale=int(yscale)

Z = image[:,(centery-yscale):(centery+yscale+1),centerx]
ZT=Z.T

ax = plt.subplot(221)
im = plt.imshow(ZT, cmap=cm.spectral,
                origin='lower', aspect='equal',
                extent=[min(z),max(z),min(y),max(y)])

xlabel('|')

Z = image[:,(centery-yscale):(centery+yscale+1),centerx]
for i in range(-xscale,xscale+1):
   for j in range(0,nz):
       Z[j,i+xscale] = image[j,centery+i,centerx+i]
ZT=Z.T
ax = plt.subplot(222)
im = plt.imshow(ZT, cmap=cm.spectral,
                origin='lower', aspect='equal',
                extent=[min(z),max(z),min(y),max(y)])
xlabel('/')

Z = image[:,centery,(centerx-xscale):(centerx+xscale+1)]
ZT=Z.T
ax = plt.subplot(223)
im = plt.imshow(ZT, cmap=cm.spectral,
                origin='lower', aspect='equal',
                extent=[min(z),max(z),min(x),max(x)])
xlabel('velocity(km/s) -')
ylabel('position')

Z = image[:,centery,(centerx-xscale):(centerx+xscale+1)]
for i in range(-xscale,xscale+1):
   for j in range(0,nz):
       Z[j,i+xscale] = image[j,centery-i,centerx+i]
ZT=Z.T
ax = plt.subplot(224)
im = plt.imshow(ZT, cmap=cm.spectral,
                origin='lower', aspect='equal',
                extent=[min(z),max(z),min(x),max(x)])
xlabel('\')

plt.show()
=============================================================

pvdiagram2.0.py

转载本文请联系原作者获取授权，同时请注明本文来自钱磊科学网博客。
链接地址：https://blog.sciencenet.cn/blog-117333-574518.html

上一篇：巨冰天外来？
下一篇：python数据处理笔记（三）通道积分图