Python Notes

The following are a few notes on the syntax, commands, etc. for python programming Assumes python version 3.x, which differs substantially from the 2.x versions. Python is not backwards-compatible, so 2.x codes don't necessarily work in 3.x.

python                   # Run python, where is a python program
python arg1 arg2         # Run python program with arguments
                                 #  inside the code you must 'import sys', and then
                                 #    str(sys.argv[0]) is ''
                                 #    str(sys.argv[1]) is the first argument, etc.
print("It is possible to break long lines of python code",  # breaking long lines
 "by preceding the carriage return with a (, \, or \

---------------Numbers, Statistics, Vectors--------------------

# + - / * do what you'd think. ** is a power, abs is absolute value
# pow(2,7) is the same as 2**7, max and min do what you'd think
# round(2.71828,4) gives 2.7183
# 10%3         # 10mod(3) = 1
# str(10%3)    # '1' returned as a string
# float('3.44') # converts a string to a floating point number
                # NOTE: ALL PYTHON inputs from the user are interpreted as strings

# Running it from the command line; most useful as a calculator
python                           # Start python (or python3 for version 3.3)
>>> import numpy as np           # import package numpy
>>> print(np.sin(np.pi/4.))      # prints 0.7071067...; need np.sin as sin is not defined
>>> a=500.3
>>> b=3
>>> a/b                          # prints out 166.7666
>>> _/5.                         # prints 33.3533   _ is the result from the last command
>>> Ctl-D                        # quit

x=[1.3,2.4,8]           # is a list, these behave like strings
x=(1.3,2.4,8)           # is a tuple, these are like lists that can't be changed
x=np.array([1.3,2.4,8]) # is what you want for a vector (1D array)
x.T                     # transpose of x, where x is a numpy array
np.zeros(10,  # initialize an integer array of 10 elements to zero
np.ceil([2.87])   # rounds up to 3
np.array([2.87,4.8], dtype=int) # rounds to [2,4]  !!Use this, and not!!
np.floor([2.87,4.8]) # rounds to [2.,4.] but these are reals
np.arange(0,5,0.1) # makes [0, 0.1, 0.2 ... 4.9, 5.0] vector
np.linspace(0,100,5) # makes [0, 25, 50, 75, 100] vector
np.random.rand(5) # five random numbers between 0 and 1
np.random.randn(5) # five random numbers from std normal distribution
np.random.randint(0,99,5) # five random integers between 0 and 99
np.mean(x)       # mean of x
np,median(x)     # median of x
np.var(x)        # variance of x   DONT FORGET THE 'np.' 
np.std(x)        # stddev of x
np.log(x)        # ln(x)
np.log10(x)      # log10(x)
np.exp(x)        # exp(x)
np.arctan(x)     # atan(x)
np.nanmedian(x)  # Take the median of values that are not nan.
np.percentile(a,85) # value of the 85th percentile of array a. Very useful for finding 'sigma'
                    # in the presence of outliers! e.g. sigma=np.percentile(a,83.3)-np.percentile,50)

import time
t0=time.time()  #Current time. Place at beginning and end to time a section of code

import matplotlib.pyplot as plt  #MATLAB-like plots
from matplotlib import pyplot as plt # Equivalent syntax for importing
import matplotlib.ticker as ticker  #for full control over tick mark placement
plt.rcParams['font.size']=14     # make the default labels bigger

fig=plt.figure(figsize=(8,15))  #fig is a figure object, change its default shape
ax1=fig.add_subplot(122) #1 row 2 columns, access plot 2  ;  ax1 is a plotting object, a container within fig
fig, [[ax1, ax2], [ax3, ax4]] = plt.subplots(2, 2, figsize=(6, 6)) # another way to set up 4 plotting objects in a 2x2
ax1.plot(data['B-V'],data['mag'], 'bx')  #plot blue x's. r-- is red dashed lines, etc.
ax1.plot(x,y,lw=2.2)  #Thicker line width
ax1.set_ylim(-2,2) #another way to set y limits
                labelrotation=30,width=3,labelsize=10,pad=6,labelcolor='k') #tick mark size, color, angle, location
ax1.xaxis.set_major_locator(ticker.IndexLocator(base=3, offset=2)) # major ticks at 2,5,8,11, etc. other options exist
ax1.set_ylabel('Mag', fontsize=14)  # y-label
ax1.set_xlabel('B-V', fontsize=14)  # x-label
ax1.set_yscale('log')  # log_10 scale on y-axis
ax1.set_title(r'$\sigma_{15}$ Some title')  # A plot title. The r means a raw string (for latex \)
ax1.text(60, .025, r'$\mu=100,\ \sigma=15$')
ax1.hist(x,30,facecolor='b',range=[0.9,1.1]) #histogram of x, 30 bins in [0.9,1.1]
ax1.grid(True)  #add a grid
ax1.scatter(x,y,c=z,s=sz,vmin=0.8,vmax=1.2) #scatter plot, with colors defined by z, size by sz, over range 0.8-1.2

plt.subplot(325)  #3 rows, 2 columns of plots; access plot number 5 of these
plt.subplots_adjust(wspace=0.5, hspace=0.5) #extra space between the subplots
plt.subplot2grid((6, 8), (0, 3), rowspan=2, colspan=2)  # More flexible subplots. (6,8) is grid size
                                                        # (0,3) is starting location
plt.axis([-0.2,0.8,11,5])  #another way to do limits on x- and y-axes
plt.savefig('output.png') # Save to a file         # Displays the plot
fig.add_subplot(grid[-1:],sharex=ax1) # add subplot at bottom, then plt.xlabel, plt.plot etc to plot there

# colors: (r)ed, (g)reen, (b)lue, (w)hite, blac(k), (c)yan, (m)agenta, (y)ellow,
#  hundreds of others like 'darkblue', 'forestgreen', 'brown', etc.
# linestyles:  : dotted; - solid; -- dashed; -. dot-dashed;  possible to make your own
# point types: o filled circle; d or D: diamonds; s squares; + : plus ; * star ; v and ^ triangle

# Polar plot example plotting (r,theta) and using other columns for point size and labels
ax=fig.add_subplot(111,projection='polar')   # make a polar plot
 # An example using sorted data in a file to control point type
 # ax.scatter allows you to use a column to define size and color but not point type. this is a workaround
ax.set_rmax(3.2)  # radial maximum for a polar plot
ax.set_rlabel_position(313)  # Angle in degrees where the labels run
ax.set_theta_zero_location("S") # Where to put zero degrees on the plot
  # Use a label in a column to identify each point
for iv in range(np.size(data['r'])) :
  # Example of a legend
ax.legend(['1-6','6-25','26-50','>50'],fontsize=14,loc='upper right', bbox_to_anchor=(1.15, 1.0))
plt.title('Number of O Stars in Nearby Associations',fontsize=20)

import matplotlib.pyplot as plt
import numpy as np
from matplotlib.animation import FuncAnimation
# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,10))
yrs=120           #length of simulation in years
tinterval=15./365.25  #time step in years #number of frames in the animation
Re=1                           #Earth radius
Rj=5.2                         #Jupiter radius
Rs=9.8                         #Saturn radius
pjup=11.86     #Jupiter period
psat=29.46     #Saturn period
thear=(2*np.pi)*tinterval #Theta changes per interval
t0e=0.*(np.pi/180.) #Initial thetas

ax.plot([0],[0],'yo',ms=16) #Sun

th=np.arange(0,365,5)*np.pi/180. #plot orbits
ax.plot(th, rade,'b',lw=0.7)
ax.plot(th, radj,'k',lw=0.7)
ax.plot(th, rads,'k',lw=0.7)
ax.text(0, 10.25, 'C-0')

 #animated parts to plot
jupiter, = ax.plot([],[], 'o',mfc='indianred',ms=10)
saturnr, = ax.plot([],[], 'k_',ms=18) #add a ring to saturn as it goes around
saturn, = ax.plot([],[], 'ko',ms=9)
earth, = ax.plot([],[], 'bo',ms=7)
time_text = ax.text(0.02, 0.95, '', transform=ax.transAxes)
c1_text = ax.text(242.7*np.pi/180., 10.8, '')
c2_text = ax.text(125.4*np.pi/180., 10.6, '')
c3_text = ax.text(8.2*np.pi/180., 10.25, '')
c4_text = ax.text(250.9*np.pi/180., 10.8, '')
c5_text = ax.text(133.6*np.pi/180., 10.6, '')
c6_text = ax.text(16.3*np.pi/180., 10.25, '')

# initialization function with blank data
def init():
        #returns ax.plot commands
    return  jupiter, saturn, saturnr, earth, time_text, c1_text, c2_text, c3_text, c4_text, c5_text, c6_text

# animation function.  This is called sequentially
def animate(i):
    if((aa>19.7) & (aa<19.9)):
    if((aa>39.6) & (aa<39.8)):
    if((aa>59.4) & (aa<59.6)):
    if((aa>79.5) & (aa<79.7)):
    if((aa>99.2) & (aa<99.4)):

    the = i*thear+t0e
    thj = i*thjup+t0j
    ths = i*thsat+t0s

    rad = [Re,Rj,Rs]
    th = [the,thj,ths]
    earth.set_data(the, Re)
    time_text.set_text('%.2f years' % aa)
        # return ax.plot commands
    return jupiter, saturn, saturnr, earth, time_text, c1_text, c2_text, c3_text, c4_text, c5_text, c6_text

# call the animator.  blit=True means only re-draw the parts that have changed.
anim = FuncAnimation(fig, animate, init_func=init, frames=d, interval=20, blit=True)'test.mp4', fps=30, extra_args=['-vcodec', 'libx264'])

---------------Displaying and Altering Images--------------------
plt.figure(figsize=(20,10))  #set the size
plt.axis('off')       # shuts off drawing of x,y axes 
plt.imshow(imname[20:42,29:48],origin='lower',vmin=-120,  #imname is np.array, vmin and vmax are scale limits
  vmax=20,cmap='jet',interpolation='nearest',aspect='equal') # cmap is colormap, aspect=equal gives square pixels

plt.imshow(blurredz[20:42,29:48],origin='lower',vmin=-120,vmax=20, # OR try using 'extent' to fix the labels
  cmap='jet',interpolation='nearest',       # Fits image [x,y] --> np.array [y,x] via astropy
  aspect='equal',extent=(29,48,20,42))      # But imshow displays np.array [y,x] as [x,y] again (I think)
                                            # Note, labels x,y flip for extent. 
                                            # If you were to put [2900,4800,2000,4200] for extent it would only
                                            # change the labels

blurredz = gaussian_filter(data3, sigma=1.5, mode='nearest') # smooth an image with a gaussian filter


# Assigning, concantating strings
ss=c+" "+c      # ss = "hello hello" 

# uppercase and lowercase
ss= "HEllo"
ss.upper()     #would be HELLO
ss.lower()     #would be hello
ss.isupper()   #returns False because HEllo is not all uppercase

# string length
len(ss)       # equals 6

# substrings and positions in a string, replacing strings
ss=Hi there"
ss[1:4]      # equals 'i th' because the index starts at 0
ss[len(ss)-3:len(ss)-1]  # 'ere' the last three characters
ss[-2]       # second to last character, 'r'
ss[-2:]      # second to last character to the end of the string, 're'
ss[::-1]     #reverses the string, 'ereht iH'
ss.index("h")    # equals 4 in this case
ss.index("her") # also equals 4, the location of the beginning of 'her' 
ss.replace("her","eak")  # gives 'Hi teake'

rstr="10:43:24.5800 -59:36:03.300"
ram=float((rstr.split(" ",1)[0]).split(":",3)[1])  # sets ram=43.0

#         Formatted print statements
print("{0:5.3f} {1:3d}".format(var,num))  #the 0: is the first argument (var),
                                          #    1: is the second argument (num)
print('Input filename: ',end="")          # stops the carriage return
print("{0:02.0f}:{1:02.0f}:{2:07.4f}".format(rah,ram,ras)) # does leading zeros for RA, e.g. 09:08:04.233
                                              # Be careful, as this format rounds up or down

print('12'.zfill(5))             # for leading zeros, prints 00012
if(np.mod(ix, :
   print("{0:s}Doing the composite median {1:2d}%".format("\r",*ix/ixbig)),end="")
                                 # prints out progress of task in a loop, every 5% or so,
                                 # overwrites the 5% in place with \r format

-----I/O into a program, open and write to output files -------

name=input("Input filename: ")                    # read name from terminal
oname="out_"+name               # make output filename
fin= open(name,'r')             # open input filename to read  'a' would be to append
fout= open(oname,'w')           # open output filename to write (overwrite) an entirely new file
for line in fin:                # loop until run out of lines
 line=line.strip()            # get rid of invisible characters like \n
 columns=line.split()         # split the fields of the line
 var1=float(columns[0])      # columns[0] and columns[1] are strings
 fout.write('{0:9.4f} {1:9.4f}\n'.format(var1-var2,var1*var2))

tfile=open('times', "r")
line=tfile.readline()    # read in first line
while line:              # read in subsequent nonzero lines

name.readable()   # True if file can be read
name.readline()   # reads in a line
name.readlines()  # reads file into a list
name.write("\n adding this line to file") # for appending to existing files

a,b = float(input("Type in two values: ").split())  # input 2 real values

from import ascii'junk.txt', header_start=0, data_start=1, delimiter=" ")
  #junk.txt has a line of header: name, B-V, and mag, and then data
data.more()                      #To interactively page through the table                        #shows what columns are what and their type
data.colnames                    #list column names
data.rename_column('dia','D')    #rename a column
data['newcol'] = [2,4,5,2,3,5]   #add a column named 'newcol'
data.remove_column('newcol')     #remove a column 
data['dia'].format = '7.3f'      #changes format of column 'dia'
data['dia','B-V']                #a table with only these two columns
data.add_row(3344,344.4,-13.4,12.3,8.9]  #add a row
data1=data[np.where((data['long']<20) & (data['dia']>10))] #extract subset
np.where(a>4)  # returns a 'tuple', with the indices where the condition is true
               # that's why np.where is inside the [ ] in the above subset example
               # most of the time you probably don't want to be doing math on np.where
               # because you'll be doing math on the indices, not on the values
np.where(a>4, a, 5*a) # returns 'a' when a>4, and 5*a otherwise, each in a tuple
np.reshape(a,np.product(a.shape))  # converts a 2D array to a 1D array

------- FITS I/O----------
from import fits'test1.fits') # opens the HDU list
data1=hdul[0].data   # the data
hdr1=hdul[0].header   # the header
hdr1['NAXIS1']  # value of the Fits keyword NAXIS1
ixlim=data1.shape[1]  # Warning!  astropy transposes the array, so x is now y
                      # it will transpose it back again when writing out
print(data1[2:4,5:6]) #prints value of (3,6) and (4,6)

hduout=fits.PrimaryHDU(dataout.astype(np.float32)) #puts array dataout into a Fits HDU. Default type is 'double' 
hduout.writeto('test2.fits', overwrite=True) # Writes the HDU to an output file

------- if, for, while statements; looping ----------
flag2 = False
if flag1 and flag2:    # you can also use 'or'.  DONT FORGET the :
elif flag1 and not(flag2): # use of 'not' function
 print("maybe ok")
 print("not ok")
 quit()          #end program if not ok

if num1 >= num2 and num2 != num3:  # For numbers
print("done with loop")

if ((a[i,j] == 0) | (a[i,j] == 2) & (t==8)) :  # note use of ==, |, &

for i in range(0,upper) :  # 0,1,2,...upper-1  i.e., the indexed values for an array of length upper

for i in range(data.shape[0]):
  vel[i] = i*15.-150             #set some values to array vel[i]

for iv in list(range(8, 12)) :

for iv in range(np.size(vs)) :

while (i< len(data['ra'])): 
print("done with loop")

while True:
  a=int(input("Input an integer: "))
  break                              # break out of what would be an infinite while loop
 except ValueError:                  # If they didnt input an integer
  print("You need an integer here")

ixmin=np.zeros(nfiles, dtype=int)        # initialize a 1D array with zeros
np.shape(a)                              # gives the size of the array a
np.shape(a)[0]                           # the length of the first axis of the array
cgrid=np.empty((iybig,ixbig,nfiles))     # initializa an array as empty
cgrid[::]=np.NAN                         # ...and then fill it with NANs; useful for masking

list[2]         # this is "billy"
list[-1]        # this is 55.3, go from the end
list[2:]        # starts at element 2 and goes to the end:  ["billy",55.3]
list1.extend(list2) # glues lists together: ["a","b",1,2] and overwrites list1
list2.append(5)     # list2 is now [1,2,5]
list2.index(5)      # returns 2, the index location of element '5'
list2.count(5)      # number of times '5' occurs in the list
list2.insert(1,12)  # list2 is now [1,12,2,5]
list2.remove(12)    # removes first occurrence of 12
list2.pop()         # remove last element from list
list.sort()         # sorts the list
list.reverse()      # reverse sorts the list
list2.clear()       # clears the list
list3=list1.copy()  # copies a list
list4=[[0,3,5],[9,4,2],[0]]  # list4[1][0] = 9 in this case

array = (2,4)       # this is a 'tuple' and can never be changed
                    # so you can't change the 4 to a 5, for example (seems stupid)
dir(task)           # lists all the attributes of a task

dateno = {         # this defines the dictionary with a key/value
"jan": "0",        # e.g. mar/59
"feb": "31",
"mar": "59",
"apr": "90",
"may": "120",
"jun": "151",
"jul": "181",
"aug": "212",
"sep": "243",
"oct": "273",
"nov": "304",
"dec": "334",
while (True):
 a,b = input("Date (e.g. Feb 22): ").split()  #input multiple values on one line
 aa = dateno.get(a.lower()[0:3],"Bad")        # Apr -> 90, etc. Apq -> "Bad"
 if aa == "Bad":

---------------Error Handling--------------------

while True:
  a=int(input("Input an integer: "))
  break                              # break out of the infinite while loop
 except ValueError:                  # If they didnt input an integer
  print("You need an integer here")

---------------Classes --------------------
class exam:              # Define a new class of objects
   def __init__(self,question,answer):   # exam objects consist of a question and an answer = question   # will refer to question part of the exam element e1
      self.ans = answer      # e1.ans will refer to answer part of the exam element e1

                # make a list of questions
questions = [
"Is the sky blue?\n(a) yes\n(b) no\n\n",
"Does Trump have yellow hair?\n(a) yes\n(b) no\n\n",
"What color are the Rams jerseys?\n(a) purple\n(b) pink\n(c)yellow and blue\n\n",

key = [                  #the list 'key' contains exam objects

for i in key:            #each i is an exam object
 answer=input(   #asks for user input with the prompt = questions[i]
 if(answer == i.ans):
   score += 1
print("You got " + str(score) + "/" + str(len(key)) + " Correct")

---------------Subroutines (functions)--------------------
def yrs_since_HS(age):
print(yrs_since_HS(55)) # gives 37.  Input was age, returned value is age-18

---------------Modules and Importing--------------------
import numpy as np
np.pi                 # 3.14159265...
pip install module
pip uninstall module

---------------Interact with the operating system--------------------

# To execute a command from the OS and print the results to STDOUT, precede with a '!'

# To call a command from the OS:
from subprocess import call
call(["/bin/ls","-lt"])  # if successful it simply exits with code '0'

# To call a command from the OS that uses pipes:
import subprocess
p1=subprocess.Popen(["/bin/ls"],stdout=subprocess.PIPE)   #first command
p2=subprocess.Popen(["/usr/bin/awk","END{print NR}"],stdin=p1.stdout,stdout=subprocess.PIPE) #pipe to second command
p1.stdout.close() #apparently so the order of finishing doesn't matter
output=p2.communicate()[0] #set the result from the pipe (a byte string) to 'output'
print(output.decode('ascii')) #print out this byte string as ascii

# To write python commands to a file, and then execute that file in a double loop
# Note if you forget the () in g.close the python code will skip one command
#  and you get a real mess, with only 5 output files instead of 6 and mixed up contents
#  That's because the file remains open, and only closes when it reloops
from subprocess import call
import numpy as np
for iv in range(np.size(vs)) :
 for id in range(np.size(denz)) :
  fout.write('{0:10.2f}  {1:4.2f} VS DEN\n'.format(vs[iv],denz[id]))

  g.write('cp -f jj test_{0:2d}_{1:4.2f}.txt\n'.format(int(vs[iv]),denz[id]))
  g.close()   # Don't forget the ()!!!

---------------Jupyter Notebook--------------------

# To start up Jupyter:
jupyter notebook  # will create a .ipynb file (the notebook) in the current directory
                  # uses  ~/.jupyter/nbconfig/notebook.json for shortcuts
                  # These are vi-like
                  # Jupyter has an insert mode and an escape mode like vi
                  # Don't be afraid to restart the kernel if it gets confused

jupyter nbconvert --to python nbname # Converts nbname.ipynb to, a python code file
                                     # you can run via 'python'
       # Note, jupyter notebooks are really json files, meaning they
       # are text files that can be edited