In [26]:
# gdapy01_00 clears all variables and import various modules
# you must run this cell first, before running any of others
# History:
# 2022/11/05 - Created by W. Menke, using edapy_01_00 as template
# 2024/05/23 - Patches to accomodate new verions of numpy, scipy, matplotlib
# patches dot product to return scalar value
# A reset deletes any previously-created variables from memory
%reset -f
# import various packages
import numpy as np # matrices & etc
from matplotlib import pyplot as plt # general plotting
from matplotlib import patches # plotting shapes
from datetime import date # date and time arithmetic
from math import exp, pi, sin, sqrt # math functions
import scipy.linalg as la # linear algebra functions
import os # operating system
# function to make a numpy N-by-1 column vector
# c=eda_cvec(v1, v2, ...) from a list of several
# array-like entities v1, v2, including a number
# a list of numbers, a tuple of numbers, an N-by-0 np array
# and a N-by-1 np array. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it
# It promotes integers to floats, and integers
# and floats to complex, by context.
# This version concatenates many argments,
# whereas c=eda_cvec1(v1) takes just one argiment.
# I recommend always using eda_cvec(v1, v2, ...)
def gda_cvec(*argv):
t = int;
Nt = 0;
for a in argv:
v = gda_cvec1(a);
N,M = np.shape(v);
Nt = Nt + N;
if( N==0 ):
continue; # skip vector of zero length
if (t==int) and isinstance(v[0,0],float):
t=float;
elif isinstance(v[0,0],complex):
t=complex;
w = np.zeros((Nt,1),dtype=t);
Nt = 0;
for a in argv:
v = gda_cvec1(a);
N,M = np.shape(v);
w[Nt:Nt+N,0:1] = v; # patch 20230418 was #w[Nt:Nt+N,0] = v[0:N,0];
Nt = Nt + N;
return w;
# function to make a numpy N-by-1 column vector
# c=gda_cvec1(v) from entity v that is array-like,
# including a number, a list of numbers, a tuple
# of numbers, an N-by-0 np array and a N-by1 np array.
# It promotes integers to floats, and integers
# and floats to complex, by context. The function
# also insures that, if v is an np array, that
# c is a copy, as contrasted to a view, of it.
# This version takes just one input argment.
# whereas c=gda_cvec(v1,v2,...) concatenates
# many argiments.
def gda_cvec1(v):
if isinstance(v, int) or isinstance(v, np.int32):
w = np.zeros((1,1),dtype=int);
w[0,0] = v;
return w;
elif isinstance(v, float):
w = np.zeros((1,1),dtype=float);
w[0,0] = v;
return w;
elif isinstance(v, complex):
w = np.zeros((1,1),dtype=complex);
w[0,0] = v;
return w;
elif isinstance(v, np.ndarray):
s = np.shape(v);
if len(s) == 1:
return np.copy(np.reshape(v,(s[0],1)));
else:
[r,c]=s;
if( c==1 ):
return(np.copy(v));
elif(r==1):
return(np.copy(v.T));
else:
raise TypeError("gda_cvec: %d by %d ndarray not allowed" % (r, c));
elif isinstance(v, list):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: list contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(v); # patch v -> np.array(v)
return w;
elif isinstance(v, tuple):
r = len(v);
t = int;
for vi in v:
if isinstance(vi,int) or isinstance(vi, np.int32): #patch v->vi 20230418
pass;
elif isinstance(vi,float):
t=float;
elif isinstance(vi,complex):
t=complex;
break;
else:
raise TypeError("gda_cvec: tuple contains unsupported type %s" % type(vi));
w = np.zeros((r,1),dtype=t);
w[:,0] = np.array(list(v)); # patch v -> np.array(list(v));
return w;
else:
raise TypeError("gda_cvec: %s not supported" % type(v));
In [27]:
# gdapy01_01
# displays the current date
thedate = date.today();
print(thedate);
2024-05-23
In [28]:
# gdapy01_02
# examples of directory manipulation, support Section I.3
# print working directory
# get current directory, with is presumably gdapy/Notebooks
gdapy_dir=os.getcwd();
print(gdapy_dir);
# change directory to one up
os.chdir("..");
# print current working directory (folder)
thisdir=os.getcwd();
print(thisdir);
# change to Notebooks folder
os.chdir('Notebooks');
thisdir=os.getcwd();
print(thisdir);
# print directory listing of notebook folder
mydirlist=os.listdir();
print(mydirlist);
# change back to original folder
os.chdir(gdapy_dir);
c:\bill\gda5ed\PatchPyton20240523\gdapy\Notebooks c:\bill\gda5ed\PatchPyton20240523\gdapy c:\bill\gda5ed\PatchPyton20240523\gdapy\Notebooks ['.ipynb_checkpoints', '.virtual_documents', 'gdapy_01.ipynb', 'gdapy_02.ipynb', 'gdapy_03.ipynb', 'gdapy_04.ipynb', 'gdapy_05.ipynb', 'gdapy_06.ipynb', 'gdapy_07.ipynb', 'gdapy_08.ipynb', 'gdapy_09.ipynb', 'gdapy_10.ipynb', 'gdapy_11.ipynb', 'gdapy_12.ipynb', 'gdapy_13.ipynb', 'gdapy_14.ipynb', 'gdapy_15.ipynb', 'gdapy_16.ipynb', 'gdapy_17.ipynb', 'junk.ipynb']
In [29]:
# gdapy01_03: simple algebra
a = 4.5;
b = 5.1;
c = a+b;
print(c);
9.6
In [30]:
# gdapy01_04: algebra with square root and power
a = 6.0;
b = 8.0;
c = sqrt(a**2 + b**2);
print(c);
10.0
In [31]:
# gdapy01_05: algebra with sine
n = 3.0; x = 4.0; x0 = 1.0; L = 6.0;
c = sin(n*pi*(x-x0)/L);
print(c);
-1.0
In [32]:
# gdapy01_06 list, tupples, vectors and matrices
print("List");
mylist = [4.0, 7.0, 2.0, 6.0];
b = mylist[2];
N=len(mylist);
print("mylist", mylist);
print("b=mylist[3]", b);
print("len(mylist)",N);
print("Tuple");
mytuple = (4.0, 7.0, 2.0, 6.0);
b = mytuple[2];
N=len(mytuple);
print("mytuple", mytuple);
print("mytuple[0]", mytuple[0]);
print("len(mytuple)",N);
# create row and column vectors
# note that decimal points are needed,
# else you will get an integer vector
# which is usually not what you want
r = np.array([[2.0, 4.0, 6.0]]);
c = np.array([[1.0], [3.0], [5.0]]);
print("row-vector r");
print(r);
print("column-vector c");
print(c);
# size of vector
N,M = np.shape(r);
print("r is", N, " by ", M);
N,M = np.shape(c);
print("c is", N, " by ", M);
# element access
a=r[0,1];
b=c[2,0];
print("a=r[0,1]");
print(a);
print("b=c[2,0]");
print(b);
# example of transposition
rT = r.T;
print("rT is a column-vector");
print(rT);
cT = c.T;
print("cT is a row-vector");
print(cT);
# vectors of zeros
N=3;
r = np.zeros((1,N));
c = np.zeros((N,1));
print("row-vector r of zeros");
print(r);
print("column-vector c of zeros");
print(c);
# vectors of ones
N=3;
r = np.ones((1,N));
c = np.ones((N,1));
print("row-vector r of ones");
print(r);
print("column-vector c of ones");
print(c);
# convert column-vector into Nx0 vector v
v = c.ravel();
print("Nx0 vector");
print(v);
# gda_cvec() function
c = gda_cvec(2.0, 4.0, 6.0);
print("column-vector c");
print(c);
# square matrix A
# note that decimal points are needed,
# else you will get an integer vector
# which is usually not what you want
A = np.array( [ [1., 4., 7.], [2., 5., 8.], [3., 6., 9.]] );
print("A");
print(A);
# find shape of A
N, M = np.shape(A);
print('shape of A:',N,M);
# 3x3 square matrix A of zeros
N=3;
M=3;
B = np.zeros((N,M));
print("B");
print(B);
# 3x3 square matrix B of ones
N=3;
M=3;
C = np.ones((N,M));
print("C");
print(C);
# element acces
s = A[0,2];
print("s = A[0,2]");
print(s);
# define matrix A
A = np.array( [ [1., 0., 2.], [0., 1., 0.], [2., 0., 1.]] );
print("A");
print(A);
# define matrix B
B = np.array( [ [1., 0., -1.], [0., 2., 0.], [1., 0., 3.]] );
print("B");
print(B);
S = A+B;
print("S = A+B");
print(S);
D = A-B;
print("D = A-B");
print(D);
List mylist [4.0, 7.0, 2.0, 6.0] b=mylist[3] 2.0 len(mylist) 4 Tuple mytuple (4.0, 7.0, 2.0, 6.0) mytuple[0] 4.0 len(mytuple) 4 row-vector r [[2. 4. 6.]] column-vector c [[1.] [3.] [5.]] r is 1 by 3 c is 3 by 1 a=r[0,1] 4.0 b=c[2,0] 5.0 rT is a column-vector [[2.] [4.] [6.]] cT is a row-vector [[1. 3. 5.]] row-vector r of zeros [[0. 0. 0.]] column-vector c of zeros [[0.] [0.] [0.]] row-vector r of ones [[1. 1. 1.]] column-vector c of ones [[1.] [1.] [1.]] Nx0 vector [1. 1. 1.] column-vector c [[2.] [4.] [6.]] A [[1. 4. 7.] [2. 5. 8.] [3. 6. 9.]] shape of A: 3 3 B [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.]] C [[1. 1. 1.] [1. 1. 1.] [1. 1. 1.]] s = A[0,2] 7.0 A [[1. 0. 2.] [0. 1. 0.] [2. 0. 1.]] B [[ 1. 0. -1.] [ 0. 2. 0.] [ 1. 0. 3.]] S = A+B [[2. 0. 1.] [0. 3. 0.] [3. 0. 4.]] D = A-B [[ 0. 0. 3.] [ 0. -1. 0.] [ 1. 0. -2.]]
In [33]:
# gdapy01_07
# examples of vector and matrix algebra
# column vector c
a = gda_cvec(1.,3.,5.);
print("a");
print(a);
# column vector r
b = gda_cvec(2.,4.,6.);
print("b");
print(b);
# square matrix A
A = np.array( [[1., 0., 2.],
[0., 1., 0.],
[2., 0., 1.]] );
print("A");
print(A);
# square matrix B
B = np.array( [[1., 0., -1.],
[0., 2., 0.],
[-1., 0., 3.]] );
print("B");
print(B);
# identity matrix
I = np.identity(3);
print('I');
print(I);
# inner product: s = a-transpose time b
s = np.matmul(a.T,b); s=s[0,0];
print("s");
print(s);
# outer product: T = a times b-transpose
T = np.matmul(a,b.T);
print("T");
print(T);
# matrix times a vector: P = A times B
c = np.matmul(A,b);
print("c=A*b");
print(c);
# matrix times a matrix: P = A times B
P = np.matmul(A,B);
print("P");
print(P);
# d = element-wise multiply of a and b
d = np.multiply(a,b);
print("d");
print(d);
a [[1.] [3.] [5.]] b [[2.] [4.] [6.]] A [[1. 0. 2.] [0. 1. 0.] [2. 0. 1.]] B [[ 1. 0. -1.] [ 0. 2. 0.] [-1. 0. 3.]] I [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] s 44.0 T [[ 2. 4. 6.] [ 6. 12. 18.] [10. 20. 30.]] c=A*b [[14.] [ 4.] [10.]] P [[-1. 0. 5.] [ 0. 2. 0.] [ 1. 0. 1.]] d [[ 2.] [12.] [30.]]
In [34]:
# gdapy01_08
# sub-matrices using colon operator
# column-vector a
a = gda_cvec([1.,2.,3.]);
print("a");
print(a);
# top 2 elements of vector
x = a[0:2,0:1];
print("x");
print(x);
# square matrix M
B = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("B");
print(B);
# column of matrix
y=B[0:3,1:2];
print("y is 2nd col of B");
print(c);
# row of matrix
z=B[1:2,0:3];
print("d is 2nd row of B");
print(d);
# block of a matrix
T=B[1:3,1:3];
print("T is lower-right 2x2 block of B");
print(T);
# if you want them independent, you must copy them
x = np.copy( a[0:2,0:1] );
y = np.copy( B[0:3,1:2] );
c = np.copy( B[1:2,0:3] );
T = np.copy( B[1:3,1:3] );
# make a column-vector
c = gda_cvec( 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.);
print("column-vector c");
print(c);
d = c[1:10:2,0:1];
print("d is subsampled");
print(d);
e = c[10:5:-1,0:1];
print("e is flipped");
print(e);
a [[1.] [2.] [3.]] x [[1.] [2.]] B [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] y is 2nd col of B [[14.] [ 4.] [10.]] d is 2nd row of B [[ 2.] [12.] [30.]] T is lower-right 2x2 block of B [[5. 6.] [8. 9.]] column-vector c [[ 0.] [ 1.] [ 2.] [ 3.] [ 4.] [ 5.] [ 6.] [ 7.] [ 8.] [ 9.] [10.] [11.] [12.]] d is subsampled [[1.] [3.] [5.] [7.] [9.]] e is flipped [[10.] [ 9.] [ 8.] [ 7.] [ 6.]]
In [35]:
# edapy_01_09
# t-axis
N=6;
Dt = 2.0;
t = gda_cvec(np.linspace(0.0,Dt*(N-1),N));
print("time axis t");
print(t);
time axis t [[ 0.] [ 2.] [ 4.] [ 6.] [ 8.] [10.]]
In [36]:
# gdama01_10
# matrix inverse, determinant, eigenvalues
# square matrix A
A = np.array( [[1., 2., 3.],
[5., 7., 11.],
[13., 17., 19.]] );
print("A");
print(A);
# vector b
b = gda_cvec(1.0,2.0,3.0);
print("b");
print(b);
# determinant of amatric
d = la.det(A);
print("d = det(A)");
print(d);
# inverse of matrix
B = la.inv(A);
print("B=inv(A)");
print(B);
# check inverse
C = np.matmul(A,B);
print("C = A*inv(A)");
print(C);
# check inverse
D = np.matmul(B,A);
print("D = inv(A)*A");
print(D);
# solve b = A*c via solve()
print(A);
print(b);
c = la.solve(A,b);
print("c is solution to Ac=b");
print(c);
# error
e = b - np.matmul(A,c);
print("error e=B-Ac");
print(e);
# square matrix B
B = np.array( [[1., 2., 0.],
[3., 3., 0.],
[4., 2., 4.]] );
print("B");
print(B);
# solve B=AD via solve()
D=la.solve(A,B);
print("D is solution to B=AD");
print(D);
# check result
E = B-np.matmul(A,D);
print("error E=B-A*D");
print(E);
# solve B=DA via solve()
D=la.solve(A.T,B.T).T;
print("D is solution to B=DA");
print(D);
# check result
E = B-np.matmul(D,A);
print("error E=B-D*A");
print(E);
# square symetric matrix A
A = np.array( [[1., 2., 0.],
[2., 2., 0.],
[0., 0., 4.]] );
print("A");
print(A);
# eigenvalue decomposition. Note: use eigh() not eig().
LAMBDA,V = la.eigh(A);
print("V");
print(V);
print("LAMBDA");
print(LAMBDA);
# check orthigonality
VVT = np.matmul(V,V.T);
print("VVT");
print(VVT);
# check orthogonality
VTV = np.matmul(V.T,V);
print("VTV");
print(VTV);
# reconstruct matrix
B = np.real(np.matmul( V, np.matmul( np.diag(LAMBDA), V.T)));
print("B = V*LAMBDA*VT");
print(B);
A [[ 1. 2. 3.] [ 5. 7. 11.] [13. 17. 19.]] b [[1.] [2.] [3.]] d = det(A) 23.999999999999993 B=inv(A) [[-2.25 0.54166667 0.04166667] [ 2. -0.83333333 0.16666667] [-0.25 0.375 -0.125 ]] C = A*inv(A) [[ 1.00000000e+00 1.11022302e-16 -2.77555756e-17] [-8.04911693e-16 1.00000000e+00 3.05311332e-16] [-2.35922393e-15 5.55111512e-16 1.00000000e+00]] D = inv(A)*A [[ 1.00000000e+00 2.49800181e-16 -1.38777878e-16] [ 2.77555756e-16 1.00000000e+00 6.10622664e-16] [ 1.38777878e-16 2.77555756e-17 1.00000000e+00]] [[ 1. 2. 3.] [ 5. 7. 11.] [13. 17. 19.]] [[1.] [2.] [3.]] c is solution to Ac=b [[-1.04166667] [ 0.83333333] [ 0.125 ]] error e=B-Ac [[ 0.00000000e+00] [-8.88178420e-16] [-1.77635684e-15]] B [[1. 2. 0.] [3. 3. 0.] [4. 2. 4.]] D is solution to B=AD [[-0.45833333 -2.79166667 0.16666667] [ 0.16666667 1.83333333 0.66666667] [ 0.375 0.375 -0.5 ]] error E=B-A*D [[ 0.00000000e+00 0.00000000e+00 2.77555756e-16] [ 0.00000000e+00 -2.66453526e-15 2.77555756e-16] [ 0.00000000e+00 -8.88178420e-16 8.88178420e-16]] D is solution to B=DA [[ 1.75000000e+00 -1.12500000e+00 3.75000000e-01] [-7.50000000e-01 -8.75000000e-01 6.25000000e-01] [-6.00000000e+00 2.00000000e+00 -1.11022302e-16]] error E=B-D*A [[-6.66133815e-16 -8.88178420e-16 -1.05471187e-15] [ 0.00000000e+00 0.00000000e+00 0.00000000e+00] [ 4.44089210e-16 2.22044605e-16 2.22044605e-15]] A [[1. 2. 0.] [2. 2. 0.] [0. 0. 4.]] V [[ 0.78820544 0.61541221 0. ] [-0.61541221 0.78820544 0. ] [ 0. 0. 1. ]] LAMBDA [-0.56155281 3.56155281 4. ] VVT [[1.00000000e+00 2.60380353e-17 0.00000000e+00] [2.60380353e-17 1.00000000e+00 0.00000000e+00] [0.00000000e+00 0.00000000e+00 1.00000000e+00]] VTV [[ 1.00000000e+00 -2.60380353e-17 0.00000000e+00] [-2.60380353e-17 1.00000000e+00 0.00000000e+00] [ 0.00000000e+00 0.00000000e+00 1.00000000e+00]] B = V*LAMBDA*VT [[1. 2. 0.] [2. 2. 0.] [0. 0. 4.]]
In [37]:
# gdapy01_11
# character string formatting amd lists of character strings
myfilename="mydata.txt";
print(myfilename);
# format an integer
i=1;
myfilename="myfile_%d.txt" % (i);
print(myfilename);
# format a fractioanl number
x=10.40;
mysentence = "position x is %.2f meters" % (x);
print(mysentence);
# define list f character strings
mycolorlist = ["red", "crimson", "chartruse", "teal"];
#access one element of the list
myfavoritecolor = mycolorlist[3];
print("myfavoritecolor");
print(myfavoritecolor);
mydata.txt myfile_1.txt position x is 10.40 meters myfavoritecolor teal
In [38]:
# gdapy01_12
# example of a simple for-loop used to extract the diagonal of a matrix
# square matrix A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# for loop that extracts the main diagonal of a matrix
a = np.zeros((3,1));
for i in range(3):
a[i,0]=A[i,i];
print("a\n",a)
A [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] a [[1.] [5.] [9.]]
In [39]:
# gdapy01_13
# example of two nested for loops used to rearrange (flip left-right)
# the elements of a matrix
# square matrix A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A\n",A);
# create an square matrix B of zeros
B = np.zeros((3,3));
# use for loop to rearrange A into B
for i in range(3): # loop over all rows
for j in range(3): # loop over all cols
B[i,2-j]=A[i,j]; # left-right-flip
print("B\n",B);
A [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] B [[3. 2. 1.] [6. 5. 4.] [9. 8. 7.]]
In [40]:
# gdapy01_14
# example of a faor loop contain a conditional (if) statment
# create a vector
a = gda_cvec( [1.0, 2.0, 1.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );
N, M = np.shape(a);
# for loop that "clips" a vector; that is, sets elements large than 6 to 6
b = np.zeros((N,1));
for i in range(N):
if a[i,0] >= 6.:
b[i,0]=6.;
else:
b[i,0]=a[i,0];
# side-by-side concatenation vectors a and b for aeshetic printing
c = np.concatenate((a,b),axis=1);
print("a and b\n",c);
a and b [[1. 1.] [2. 2.] [1. 1.] [4. 4.] [3. 3.] [2. 2.] [6. 6.] [4. 4.] [9. 6.] [2. 2.] [1. 1.] [4. 4.]]
In [41]:
# gdapy_01_15: avoiding for-loops
# square matric A
A = np.array( [[1., 2., 3.],
[4., 5., 6.],
[7., 8., 9.]] );
print("A");
print(A);
# extract diagonal of matcix into a vector
a = gda_cvec( np.diagonal(A) );
print("a");
print(a);
# B is A flipped left to right
B = np.fliplr(A);
print("B");
print(B);
# column vector a of lenth 12
a = gda_cvec( [1.0, 2.0, 10.0, 4.0, 3.0, 2.0, 6.0, 4.0, 9.0, 2.0, 1.0, 4.0] );
print("a");
print(a)
# alternative method #1 of cliping a vector
b=np.copy(a); # start with c that is a copy of a
b[a>=6.]=6.;
print("b");
print(b);
# alternative method #1 of cliping a vector
b = np.multiply(a, a<6.) + np.multiply(6., a>=6.);
print("b\n",b);
# alternative method #3 f cliping a vector
# we use the np.where() command to finds the
# rows of a that satisfy the condition
d=np.copy(a); # start with a d that is a copy of a
r = np.where(a>=6.); # r is a data structure descibed below
k = r[0]; # k is a Nx0 list of indices
d[k,0:1]=6.; # reset the rows of d that violate the condition
# explanation of where() method
# r is a tupple
print("type r\n",type(r));
print("shape r\n",np.shape(r));
# k=r[0] gives the rows satisfying the a>=6 condition. It is a 3-by-0 np.array
print("type k=r[0]\n",type(k));
print("shape k\n",np.shape(k));
print("k\n",k);
print(" ");
print("d\n",d);
A [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] a [[1.] [5.] [9.]] B [[3. 2. 1.] [6. 5. 4.] [9. 8. 7.]] a [[ 1.] [ 2.] [10.] [ 4.] [ 3.] [ 2.] [ 6.] [ 4.] [ 9.] [ 2.] [ 1.] [ 4.]] b [[1.] [2.] [6.] [4.] [3.] [2.] [6.] [4.] [6.] [2.] [1.] [4.]] b [[1.] [2.] [6.] [4.] [3.] [2.] [6.] [4.] [6.] [2.] [1.] [4.]] type r <class 'tuple'> shape r (2, 3) type k=r[0] <class 'numpy.ndarray'> shape k (3,) k [2 6 8] d [[1.] [2.] [6.] [4.] [3.] [2.] [6.] [4.] [6.] [2.] [1.] [4.]]
In [42]:
# gdapy00_16
# load and plot 1965-2020 global temperature data
# data from:
# Hansen, J., Mki. Sato, R. Ruedy, K. Lo, D.W. Lea, and M. Medina-Elizade,
# 2006: Global temperature change. Proc. Natl. Acad. Sci., 103, 14288-14293,
# doi:10.1073/pnas.0606291103.
# load rectangular, tab-separated table of numbers into a N by M matrix D
D = np.genfromtxt("../Data/global_temp.txt", delimiter="\t");
N, M = np.shape(D);
t = np.copy( D[0:N,0:1] ); # time in column 1
d = np.copy( D[0:N,1:2] ); # temperature in column 2
# plot the data
fig1 = plt.figure(1,figsize=(8,4)); # smallish figure
ax1 = plt.subplot(1, 1, 1); # only one subplot
plt.plot(t,d,"k-"); # plot d(t) with black line
plt.plot(t,d,"ro"); # plot d(t) with black line
plt.xlabel("time (years)"); # label time axis
plt.ylabel("tempeerature (deg C)"); # label data axis
plt.show(); # display plot
In [43]:
# gdapy01_17: load text file and mmke nicer plot
# load and plot 1965-2020 global temperature data
# data from:
# Hansen, J., Mki. Sato, R. Ruedy, K. Lo, D.W. Lea, and M. Medina-Elizade,
# 2006: Global temperature change. Proc. Natl. Acad. Sci., 103, 14288-14293,
# doi:10.1073/pnas.0606291103.
# load rectangular, tab-separated table of numbers into a N by M matrix D
D = np.genfromtxt("../Data/global_temp.txt", delimiter="\t");
N, M = np.shape(D);
t = np.copy( D[0:N,0:1] ); # time in column 1
d = np.copy( D[0:N,1:2] ); # temperature in column 2
# nicer plot
fig1 = plt.figure(1,figsize=(8,4)); # smallish figure
ax1 = plt.subplot(1, 1, 1); # only one subplot
plt.axis( [1960, 2025, -0.5, 1.5]); # manually set axis limits
plt.plot(t,d,"k-"); # plot d(t) with black line
plt.plot(t,d,"ro"); # plot d(t) with red circles
plt.xlabel("time (years)"); # label time axis
plt.ylabel("temperature (deg C)"); # label data axis
plt.show(); # display p
In [44]:
# gdapy01_18:
# write data to a file
# in this example, global temperature is converted to Fahrenheit
# degC to degF conversion factor
cf = 1.8;
# data from:
# Hansen, J., Mki. Sato, R. Ruedy, K. Lo, D.W. Lea, and M. Medina-Elizade,
# 2006: Global temperature change. Proc. Natl. Acad. Sci., 103, 14288-14293,
# doi:10.1073/pnas.0606291103.
# load rectangular, tab-separated table of numbers into a N by M matrix D
D = np.genfromtxt("../Data/global_temp.txt", delimiter="\t");
N, M = np.shape(D);
t = np.copy( D[0:N,0:1] ); # time in column 1
d = np.copy( D[0:N,1:2] ); # temperature in column 2
dF = cf*d;
# load the neuse.txt data set
# convert discharge from cfs to m3/s
# and write to a file neuse_metric.txt
# concatenate time, discharge into a matrix C
DF = np.concatenate( (t,dF),axis=1);
# save matrix to file of tab-separated values
np.savetxt("../TestFolder/global_temp_in_F.txt",DF,delimiter="\t");
# plot data to verify conversion
# nicer plot
fig1 = plt.figure(1,figsize=(8,4)); # smallish figure
ax1 = plt.subplot(1, 1, 1); # only one subplot
plt.axis( [1960, 2020, -1, 3]); # manually set axis limits
plt.plot(t,dF,"k-"); # plot d(t) with black line
plt.plot(t,dF,"ro"); # plot d(t) with red circles
plt.xlabel("time (years)"); # label time axis
plt.ylabel("temperature (deg C)"); # label data axis
plt.show(); # display p
In [ ]: