代码相关（python）

一个月只能修改一次次

已于 2024-06-29 22:13:29 修改

阅读量1.6k

点赞数 3

分类专栏：代码文章标签： python

于 2020-11-12 19:45:56 首次发布

本文链接：https://blog.csdn.net/qq_40226557/article/details/109524798

版权

代码专栏收录该内容

2 篇文章 1 订阅

订阅专栏

python

程序崩溃提示符
用python的时候的各个tips
Python & 天文包

程序崩溃提示符

需要自取

         *              *
        %$@            !@#
       &%*&%          ^$#!@
      <:{!+@"|&^*&%$#!@%"}#)
     #!@%$**************&%$@~
    &++&@%***  code  ***%&%*&%
    @%*#!&**   begin  **"?)^&%
     %*#!&**   now    **?%*@#
      ++&@**************@~#^
          #!*&&^$#!@%*()
            $%(_+(%@?^
         )!(#*&++&@%%!*_#+
      !@#$)(*&^$#!@%*()*&%$@
    (*&^$#!@%*()*&%$@?#^%*(*&*
  #$)(*&^$#<{}"?)^&%$@?#^%*(*&*)
 !@#$)(*&^$#!@%*&%*&%$@~#^%*(*&*)
  )(*&^$#!@%*$#*&^$#!@%*()*&%$@~
   $)(*&^$#!@%*!@*&%$@|#^%&%$@~
    (*&^$#!@%*@|_>"}#)$(*()*&%          **   **
     (*&^$#!@%*#!*&&^$#!@%*()       **          **
      )*&%$@|%*@%*&&^$#!@%*(      **      **     **
       @|#^%&%$>:}?%*@#!@%*     **       **      **
        @%*@#!@%*!@%*$#*&@*     **         **   **
        !@%*()*&:>%#!~>:"& ** **
        <:{!+@"|&^*&%$#!@%   *
       *$#!+@"*&  &%$@~<{}"

用python的时候的各个tips

矩阵

python判断某个矩阵是否满足要求

在程序中我经常要判断某个矩阵是否满足某个条件。
1.矩阵中只要有大于等于一个元素满足条件，就返回True。

import numpy as np
data = np.array([[1,2,3],[4,5,6],[7,8,9]])
if (data>3).any(): #只要data中有任何一个或多个元素大于3，就返回True
	print('True')

返回True。
2.矩阵中所有元素满足条件，才返回True。

import numpy as np
data = np.array([[1,2,3],[4,5,6],[7,8,9]])
if (data>3).all(): #只有data中所有元素都大于3，才返回True，否则返回False
	print('True')
else:
	print('False')

返回False。

python生成二维随机数

import numpy as np
mean=5#均值
sigma=1#误差
np.random.nomal(mean,sigma)#单个随机数
np.random.normal(mean,sigma,[10,10])#10x10的随机数

文件/档

python检查某个文件存不存在

import os
if os.path.exists('/Users/mmm/file.txt'):
	print('file exists')

python添加有特定字段的文件到列表

import os
items = os.listdir('/path')
newlist = []
for name in items:
	if name.endswith('.txt'): # 添加有特定拓展名的文件到newlist
		newlist.append(name)

import os
items = os.listdir('/path')
newlist = []
for word in items:
	if 'any word' in word:
		newlist.append(word) # 添加有特定字段的文件名的文件到newlist

python矩阵保存为txt文档

np.savetxt('name.txt',array,header='xxx',fmt='%e')

header：在文档最开头的注释，必须是str格式
fmt：数据格式
如果要保存一位数组为一行，用下面的方式，否则输出的文件与预期不符：

np.savetxt('name.txt',array.reshape(1,array.shape[0]),fmt='%e')

比如：array=[1,2,3,4]，不进行reshape输出为

进行reshape后输出为：

1 2 3 4

python按行读文档

f = open('xxx.cat','r')
data_in_line = f.readlines()

data_in_line是一个列表，每个元素存储了文档中的一行

python写文档

f = open('xxx.cat','w')
lines = f.write('.................')

w表示write，r表示read。如果xxx.cat之前有内容，w模式下会覆盖之前的内容。
修改某个文档的某一行或者某几行：

f = open('xxx.cat','r')
lines = f.readlines()
lines[5] = 'modify sixth line'
lines[19] = 'modify twentieth line'
f.close()
f = open('xxx.cat','w')
f.writelines(lines)
f.close()

python文档操作

# 写文档，覆盖之前文件的内容
file=open('path/name','w')
file.write('blablablabla')
# 读文档
file=open('path/name','r')
# 写文档，不覆盖之前文件的内容，从文档末尾开始写
file=open('path/name','a')
file.write('blablablabla\n'
		   'blablablabla') # \n是换行符

字符串

python用split来拆分字符串

>x = 'asdgsdg'
>y = x.split('d')

输出结果为：

y[0] = 'as'
y[1] = 'gs'
y[2] = 'g'

如果一个cat文件有很多列，比如cat文件内容为：
阿里第三节国内高傲冷酷地说
i沃尔那个女的明白你阿松i感觉呢卡拉
哦i怕好的结果可让vi i哦为哈感觉那么
我恶化那个尽可能的到i结果仍能密码

cat = open(shangbian.cat)
cat_in_line = cat.readlines()
del cat[0] # 删掉第一行
for line in cat_in_line:
	first_column = line.split()[0] # 获得第一列
	second_column = line.split()[1] # 获得第二列
	print(first_column)
	print(second_column)

结果为：

i沃尔那个女的明白你
阿松i感觉呢卡拉
哦i怕好的结果可让vi
i哦为哈感觉那么
我恶化那个尽可能的
到i结果仍能密码

python搜索字符串某个字符的位置

#string.index('x',begin,end)返回字符串从begin到end位置之间中第一次出现x的索引。
#string.rindex('x',begin,end)返回字符串从begin到end位置之间中最后一次出现x的索引。
In [1]: a='goodsn_spec_23456_2.fits'
In [2]: a.index('_')
Out[2]: 6
In [3]: a.index('_',6)
Out[3]: 6
In [4]: a.index('_',7)
Out[4]: 11
In [5]: a.rindex('_')
Out[5]: 17
In [6]: a.rindex('_',0,15)
Out[6]: 11

给字符串前/后添加字符

请参考stackoverflow
向后填充空字符

>>> a='asdf'
>>> b=a.ljust(10)
>'asdf      '

向后填充 ‘-’ 字符

>>> a='asdf'
>>> b=a.ljust(10,'-')
>'asdf------'

向前填充空字符

>>> a='asdf'
>>> b=a.rjust(10)
>'      asdf'

向前填充 ‘-’ 字符

>>> a='asdf'
>>> b=a.ljust(10,'-')
>'------asdf'

画图

python误差棒图

plt.errorbar(x,y,xerr=xerror,yerr=yerror,fmt='.',mec='r',elinewith=1,capthick=1,capsize=1)
#xerr：x轴误差；yerr：y轴误差
#fmt：数据点的样式
#mec：数据点的颜色
#elinewidth：误差棒的粗细
#capthick：误差棒上下限的长度
#capsize：误差棒上下线的粗细

python画有百分比percentage的contour轮廓图

可能很简单的东西大家就懒得写了，找了很久没找到画百分比轮廓图的程序。就是，比如画了个散点图，你想看看这些散点的百分比轮廓，密度最高的前百分之90在哪个区域内。

import numpy as np #用它的二维直方统计
import matplotlib.pyplot as plt #用它来画轮廓图
from scipy.interpolate import interp1d #用它来插出想要的位置的值

xdata = [] #你的点对应的横轴数据
ydata = [] #你的点对应的纵轴数据
xran = [] #你想用来做统计的横轴范围
yran = [] #你想用来做统计的纵轴范围
xbins = 10 #横轴你想分几个区间
ybins = 20 #纵轴你想分几个区间
h,xedges,yedges = np.histogram2d(xdata,ydata,range=np.array([xran,yran]),bins=np.array([xbins,ybins]))
#通过xran,yran,xbins,ybins你其实已经相当于把数据平面分割成了xbins*ybins个格子，
#上面的np.histogram2d就是用来统计每个格子里都有多少个数据（散点）落在里面了
#h中包含了每个格子中有多少个数据点，xedges中包含了你的横轴每个格子的坐标，yedges同理纵轴
#现在你已经把这些格子做了统计，那么怎么知道哪个范围内包含了你想要的前百分之多少的数据点呢
percent = 0.9 #取前百分之90，想要多少直接改这个值
#percent2 = 0.5 #如果你想要另一个百分比轮廓。
temp = np.sort(h.flatten())[::-1] #把h从高到低排列，取前百分之多少（如果你想取后百分之多少，画两个轮廓就行了）
count = np.sum(temp) #算算xran,yran包裹的区域内总共有多少个点
temp2 = 0 #用来判断到哪个地方就超过你想要的百分比了
for num1 in range(0,len(temp)):
    temp2 += temp[num1] #想象一座山，从最高的那一层开始统计，掐尖。从数据点最多的格子开始，一直往山下统计，直到你想要的比例
    if temp2 > percent * count:
        break
cut1,pert1 = [temp[num1-1],(temp2-temp[num1])/count] #这个格子还不到你要的比例
cut2,pert2 = [temp[num1],temp2/count] #这个格子刚超过你要的比例
#---------------------------------------------------
#temp2 = 0 #如果你想要另一个百分比轮廓，同理上一部分
#for num2 in range(0,len(temp)):
#    temp2 += temp[num2]
#    if temp2 > percent2 * count:
#        break
#cut3 = temp[num2-1]
#cut4 = temp[num2]
#---------------------------------------------------
#现在可以插值了，因为你想要前90%，这个值介于cut1和cut2之间
f1 = interp1d([pert1,pert2],[cut1,cut2],kind='linear')
#你可以选不同的插值方法，见scipy.interpolate.interp1d文档
level_want = f(percent) #这个level值，就是到时候要输入到plt.contour的levels的值
contour = plt.contour(xedges[0:xedges.shape[0]-1],yedges[0:yedges.shape[0]-1],h.T,levels=[level_want],colors=['orange','green'])
#xedges和yedges是边界，取前边界，或者后边界都行。h.T是把h转置了，也就是转了90度。因为plt.contour中的这个轮廓项是(M,N)大小的，M是纵轴的格，N是横轴的格
#https://matplotlib.org/3.5.0/api/_as_gen/matplotlib.pyplot.contour.html?highlight=contour#matplotlib.pyplot.contour
plt.clabel(contour,contour.levels,inline=True,fontsize=10,fmt={level_want:str(percent*100)+'%'},use_clabeltext=True,inline_spacing=1) #标示比例到轮廓图上

密度分布轮廓

参考stackoverflow的回答

Error

python出现错误OSError: [Errno 24] Too many open files

在终端进行下面命令：

ulimit -n # 查看可打开的最大进程数
ulimit -n 2048 # 修改系统限制的最大进程数为2048或其他

python出现Error的解决办法

除了syntax error外，类似RuntimeError, ValueError, TypeError等等很多error都属于exceptions。
运行程序出错提示为RuntimeError等exceptions时，可以用try except else finally语句来解决。
比如，今天运行代码出现了：
RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 50000.
解决如下：

try:
	statement #可能会引起RuntimeError的命令，比如可以取消下面的一行注释试试
	#raise RuntimeError
except RuntimeError: #如果try之下的命令引起了RuntimeError，那么执行except命令
	print('引起RuntimeError了，要执行except以下语句了')
	x = 1
else: #如果try之下的命令没有引起RuntimeError，那么执行else命令
	print('没有引起RuntimeError，要执行else以下语句了')
	y = 2
finally: #不管有没有引起RuntimeError，都要执行finally语句
	print('不管有没有引起RuntimeError，都要执行finally以下语句')
	z = 3

python延迟执行下一行命令

import time
time.sleep(10) #unit:second

python的传参和传值问题

Python中用"="符号给变量赋值的时候要小心
举例：

import numpy as np
a=np.array([[1,2,3],[4,5,6]])
b=a
print('a:')
print(a)
print('b:')
print(b)

给出a和b的值

a:
[[1 2 3]
 [4 5 6]]
b:
[[1 2 3]
 [4 5 6]]

接着上面的代码

b[:,1]=0
print('a:')
print(a)
print('b:')
print(b)

给出a和b的值

a:
[[1 0 3]
 [4 0 6]]
b:
[[1 0 3]
 [4 0 6]]

a和b都改变了。也就是说对于"="会改变所有曾经 $= b$ 的变量，这个传值的规则适用于矩阵(numpy.array)，列表(list)，字典({})。
那么如果你不想改变原数组a，可以用numpy中的copy函数：

import numpy as np
a=np.array([[1,2,3],[4,5,6]])
b=np.copy(a)
b[:,1]=0
print('a:')
print(a)
print('b:')
print(b)

给出a和b的值

a:
[[1 2 3]
 [4 5 6]]
b:
[[1 0 3]
 [4 0 6]]

这样就可以只改变后来的参数b，而不改变原来的参数a。
对于数值变量，例如：

a,b=1,2
c,d=a,b
c,d=5,6
print(a,b,c,d)

则不改变原变量：

(1,2,5,6)

python添加自己写的包

写一个自己的包，然后可以像numpy一样调用：

ipython
import sys
sys.path

然后把自己的包挪到这些path中的某一个就行。
然后就可以像numpy一样调用自己的包：
比如建一个包(sys.path/mycode)写了一个脚本print_of_me.py存在mycode里，内容如下：

class print_of_me(): #有没有都可以
def ppp(a):
	print(a+' of me')

然后就可以调用了

from mycode import print_of_me as pom
pom.print_of_me.ppp('ok') #pom是module，print_of_me是class，ppp是函数

给出结果

ok of me

当然没有class也没关系，只是def也可以。
然后关于__init__.py，我发现没有也没事。注意是两个杠_ _只不过连在一起了。

Python & 天文包

astropy

astropy包可以说是用python做观测天文的基础包，功能其实很强大的。一些图和表的基础用法，包括其中还有宇宙学的一些简单计算（光度距离，角距离等）。

fits 数据处理

Python处理fits的图

Python 会交换原来fits图的x、y轴
我做psf的时候出了点问题
我从SExtractor选了点源，但是根据SExtractor的X_IMAGE和Y_IMAGE去获取pixel value时，很多值都是0。但是我去查fits的ds9，ds9给出的值不是0。
后来发现astropy.io.fits.open和pyfits.open都会原fits图的交换x轴和y轴。也就是，astropy.io的x轴是dec，y轴是ra。
如何创建binary table，astropy官网有：https://docs.astropy.org/en/stable/io/fits/usage/table.html?highlight=20A#id1
https://docs.astropy.org/en/stable/io/fits/

用三个不同filter的fits图生成rgb图

(用三个二维矩阵生成RGB图,use three two-dimension images generate RGB image)
matplotlib.pyplot中的imshow可以生成rgb图，前提是需要用np.dstack把三个波段的矩阵合并成（m,n,3）的形式。
文档中说imshow的时候，输入的三位矩阵是不需要归一化的。
原文：The Normalize instance used to scale scalar data to the [0, 1] range before mapping to colors using cmap. By default, a linear scaling mapping the lowest value to 0 and the highest to 1 is used. This parameter is ignored for RGB(A) data.原文链接：https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.imshow.html?highlight=imshow#matplotlib.pyplot.imshow
但是我直接把源数据的矩阵输入imshow的时候，画出来的rgb图有问题（某个波段特别明显，某些波段几乎没有，因为某些波段的矩阵数值很大，因此很明显，小的就不明显了）于是我先把所有波段的数据归一化到[0,1]中。然后再用imshow，具体程序如下：

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
pathr = '/pathred/red.fits'
pathg = '/pathgreen/green.fits'
pathb = '/pathblue/blue.fits'
imager = fits.open(pathr)[0].data
imageg = fits.open(pathg)[0].data
imageb = fits.open(pathb)[0].data
def rescale(arr): #normalize the true value to [0,1]
    arr_min = arr.min()
    arr_max = arr.max()
    return (arr - arr_min) / (arr_max - arr_min)
imager = rescale(imager)
imageg = rescale(imageg)
imageb = rescale(imageb)
ddd = np.dstack((imager,imageg,imageb)) #combine three 2-D array into 3-D array.
plt.imshow(np.sinh(ddd1))

示例结果如下：在这里插入图片描述

Python处理图和header

因为我psfex用不了，一直卡在安装ATLAS上很多天，而我有psf的image。因此我本来打算直接把手里的psf的image套在psfex的一个输出fits中，稍微编辑一下header。然后就发现一个bug。
步骤是这样的：
1.sextractor会给一个psfex的输出的例子——default.psf。它有两层，第一层数据为空，第二层数据是一个binary table，psf的image存在PSF_MASK的field下。维数：(1,1,31,31)一个思维的矩阵。表头如下：

[1]:import numpy as np
[2]:from astropy.io import fits
[3]:psf_def=fits.open('default.psf')
[4]:psf_def[1].header
[4]:
XTENSION= 'BINTABLE'           / THIS IS A BINARY TABLE (FROM THE LDACTOOLS)    
BITPIX  =                    8 /                                                
NAXIS   =                    2 /                                                
NAXIS1  =                 3844 / BYTES PER ROW                                  
NAXIS2  =                    1 / NUMBER OF ROWS                                 
PCOUNT  =                    0 / RANDOM PARAMETER COUNT                         
GCOUNT  =                    1 / GROUP COUNT                                    
TFIELDS =                    1 / FIELDS PER ROWS                                
EXTNAME = 'PSF_DATA'           / TABLE NAME                                     
LOADED  =                   36 / Number of loaded sources                       
ACCEPTED=                   32 / Number of accepted sources                     
CHI2    =               1.4190 / Final Chi2                                     
POLNAXIS=                    0 / Number of context parameters                   
POLNGRP =                    0 / Number of context groups                       
PSF_FWHM=               2.5813 / PSF FWHM                                       
PSF_SAMP=               0.5000 / Sampling step of the PSF data                  
PSFNAXIS=                    3 / Dimensionality of the PSF data                 
PSFAXIS1=                   31 / Number of element along this axis              
PSFAXIS2=                   31 / Number of element along this axis              
PSFAXIS3=                    1 / Number of element along this axis              
TTYPE1  = 'PSF_MASK'           / Tabulated PSF data                             
TFORM1  = '961E    '                                                            
TDIM1   = '(31, 31, 1)'

2.将1中的psf第二层的数据换成我要用的psf的数据：

psf_need=fits.open('psf.fits')
psf_def[1].data['PSF_MASK'][0,0]==psf_need[0].data

这时会报错：ValueError: could not broadcast input array from shape (69,69) into shape (31,31)。
我要用的psf的image是69x69的。
3.然后开始修改defaul.psf的表头以满足69x69:

psf_def[1].header['PSFAXIS1']=69
psf_def[1].header['PSFAXIS2']=69
psf_def[1].header['TFORM1']='4761E    '
psf_def[1].header['TDIM1']='(69, 69, 1)'

然后查看一下数据的大小：

psf_def[1].data['PSF_MASK'].shape

会报错：TypeError: buffer is too small for requested array
4.它说我需要的数组的大小太大。那么我用19x19的数组试一下：（接第二步）

psf_def[1].header['PSFAXIS1']=19
psf_def[1].header['PSFAXIS2']=19
psf_def[1].header['TFORM1']='361E    '
psf_def[1].header['TDIM1']='(19, 19, 1)'

然后查看一下数据的大小：

[11]:psf_def[1].data['PSF_MASK'].shape
[11]:
(1, 1, 31, 31)

还是没有变。
5.巧就巧在，我有一次没有先赋值（第二步），直接进行了第三步（修改header），然后成功修改了数组的shape。具体如下：

[1]:import numpy as np
[2]:from astropy.io import fits
[3]:psf_def=fits.open('default.psf')
[4]:psf_def[1].header['PSFAXIS1']=19
[5]:psf_def[1].header['PSFAXIS2']=19
[6]:psf_def[1].header['TFORM1']='361E    '
[7]:psf_def[1].header['TDIM1']='(19, 19, 1)'
[8]:psf_def[1].data['PSF_MASK'].shape
[8]:
(1, 1, 19, 19)

然后，我按照上面的方法，去修改成69x69的，还是失败：

[1]:import numpy as np
[2]:from astropy.io import fits
[3]:psf_def=fits.open('default.psf')
[4]:psf_def[1].header['PSFAXIS1']=69
[5]:psf_def[1].header['PSFAXIS2']=69
[6]:psf_def[1].header['TFORM1']='4761E    '
[7]:psf_def[1].header['TDIM1']='(69, 69, 1)'
[8]:psf_def[1].data['PSF_MASK'].shape
[8]:
TypeError: buffer is too small for requested array

还是报错，也就是说这里有两个bug：1.要想改第二层数据的shape只能先改header再查看shape，如果先看shape然后再header是没有用的。2.哪怕是按照第一点所说也要求修改后的数据shape不能大于31x31。

astropy计算观测图上两点的距离(separation)

用法如下，用SkyCoord定义两个类p1和p2，然后用属性separation就可以了，注意要加单位！

from astropy.coordinates import SkyCoord
import astropy.units as u
p1 = SkyCoord(214.9060800*u.deg,52.8444033*u.deg,frame='fk5')
p2 = SkyCoord(214.90607961*u.deg,52.84440246*u.deg,frame='fk5')
sep = p1.separation(p2)
print(sep)
$ 0d00m00.00314065s
print(sep.arcsec) #用角秒来表示两点之间的距离
$ 0.0031406469964207367
print(sep.deg) #用度来表示两点之间的距离
$ 8.724019434502046e-07
print(sep.arcmin) #用角分来表示两点之间的距离
$ 5.234411660701228e-05

用astropy.io.fits生成表格

逻辑就是先生成数据的矩阵，然后定义每一列，然后再生成表格

from astropy.io import fits
import numpy as np
array = np.array([1,2,3,4,5,6) #数据
col = fits.Column(name='id',array=array,format='K') #fits的一列
hdu = fits.BinTableHDU.from_columns([col]) #用列形成表格

前面的过程是你的每一列都只有一个数据，假如你想让fits表格中的每一行的这一列是一个数组（有多个数据）而非一个数，那么只需要改format就行了（加个P）。format参数的各个字母的含义在这里。

from astropy.io import fits
import numpy as np
array = np.array([1,2,3,4,5,6) #数据
col = fits.Column(name='id',array=array,format='PK') #fits的一列
hdu = fits.BinTableHDU.from_columns([col]) #用列形成表格

宇宙学参数

具体用法参考astropy官网

>from astropy.cosmology import FlatLambdaCDM
>h0 = 70 # km/s/Mpc
>omega_matter = 0.3
>cosmo = FlatLambdaCDM(H0=h0, Om0=omega_matter)
>#计算红移1的光度距离
>z = 1
>print(cosmo.luminosity_distance(z))
6607.657611774941 Mpc
>print(cosmo.luminosity_distance(z).to_value())
6607.657611774941

photutils

photutils包可以用来图像测光等。

photutils在图像上定义圈和环

以像素为单位的圈

from photutils.aperture import CircularAperture as ca
position = [50,50] #以(50,50)这个点为中心
aperture = ca(position,r=5) #半径是5个pixel

以图像的天球坐标wcs为单位的圈

from photutils.aperture import SkyCircularAperture as sca
from astropy.coordinates import SkyCoord
import astropy.units as u
ra,dec = 214.9060800,52.8444033
position=SkyCoord(ra*u.deg,dec*u.deg,frame='fk5') #以天球上(214.9060800,52.8444033)这个点为中心
aperture=sca(position,r=0.5*u.arcsec) #半径是0.5个角秒

这个时候你要应用到图上还需要根据图像的wcs，把aperture转换成pixel信息（因为你在之后做测光的时候是对图像上的pixel做的操作，在矩阵中只认pixel不认天球坐标）

from astropy import wcs
from astropy.io import fits
img = fits.open('example.fits')[0]
w=wcs.WCS(img.header) #提取图像wcs信息
aperture = aperture.to_pixel(w) #把原aperture的天球坐标信息转换成这张图像的pixel信息