AnalisisDeDatos

August 20, 2017

In [19]: % matplotlib inline

In [2]: from netCDF4 import Dataset

import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
import datetime

In [3]: gistemp=Dataset('gistemp1200_ERSSTv4.nc','r')# Llama netCDF en modo lectura

In [22]: print gistemp.variables

OrderedDict([(u'lat', <type 'netCDF4._netCDF4.Variable'>

float32 lat(lat)
standard_name: latitude
long_name: Latitude
units: degrees_north
unlimited dimensions:
current shape = (90,)
filling off
), (u'lon', <type 'netCDF4._netCDF4.Variable'>
float32 lon(lon)
standard_name: longitude
long_name: Longitude
units: degrees_east
unlimited dimensions:
current shape = (180,)
filling off
), (u'time', <type 'netCDF4._netCDF4.Variable'>
int32 time(time)
long_name: time
units: days since 1800-01-01 00:00:00
bounds: time_bnds
unlimited dimensions:
current shape = (1650,)
filling off
), (u'time_bnds', <type 'netCDF4._netCDF4.Variable'>
int32 time_bnds(time, nv)

1
unlimited dimensions:
current shape = (1650, 2)
filling off
), (u'tempanomaly', <type 'netCDF4._netCDF4.Variable'>
int16 tempanomaly(time, lat, lon)
long_name: Surface temperature anomaly
units: K
scale_factor: 0.01
cell_methods: time: mean
_FillValue: 32767
unlimited dimensions:
current shape = (1650, 90, 180)
filling off
)])

0.1 Definicin de la variables

In [23]: tempanomaly = np.array (gistemp.variables['tempanomaly'][:])
lat = np.array (gistemp.variables['lat'][:])
lon = np.array (gistemp.variables['lon'][:])

print len(lat)
print len(lon)
print tempanomaly.shape

90
180
(1650, 90, 180)

In [24]: time = np.array (gistemp.variables['time'][:])

time = time.astype(float)

print len(time)
print time.shape
print time.dtype

1650
(1650,)
float64

In [29]: fechas = []

for i in range(len(time)):
#print i
#print time[i]
fechas.append(datetime.datetime(1800,01,01)+ \

2
 datetime.timedelta(days = time[i]))

fechas = np.array(fechas)
print fechas

[datetime.datetime(1880, 1, 15, 0, 0) datetime.datetime(1880, 2, 15, 0, 0)

datetime.datetime(1880, 3, 15, 0, 0) ...,
datetime.datetime(2017, 4, 15, 0, 0) datetime.datetime(2017, 5, 15, 0, 0)
datetime.datetime(2017, 6, 15, 0, 0)]

0.1.1 Elimina los datos faltantes en la matriz de temperatura

In [26]: tempanomaly[tempanomaly==32767]=np.nan

0.1.2 Seleccin de un solo mapa de toda la matriz

In [27]: a = np.where(fechas == datetime.datetime(1990,07,15))
print a
print a[0][0]

Mapa = tempanomaly[a[0][0],:,:]
print Mapa.shape

(array([1326]),)
1326
(90, 180)

0.1.3 Seleccin de la serie de un pixelde toda la matriz

In [28]: print lat[45]
print lon[50]
Serie=tempanomaly[:,45,50]

Serie2=tempanomaly[:,0,50]

1.0
-79.0

0.1.4 Graficacin Mapa

In [33]: fig = plt.figure(figsize=(10,8))
ax = fig.add_subplot(111)

# Basemap es el paquete que dibuja las lneas del mapa

m = Basemap(llcrnrlat=np.min(lat),urcrnrlat=np.max(lat), \
llcrnrlon=np.min(lon),urcrnrlon=np.max(lon),\

3
 rsphere=6371200.,resolution='l',area_thresh=10000)

ny = lat.shape[0]; nx = lon.shape[0]
lons, lats = m.makegrid(nx, ny)
x,y = m(lons, lats)

cs = m.contourf(x,y,Mapa,cmap='jet')

m.colorbar(location='bottom',pad="10%")

m.drawparallels(np.arange(-90.,90,30.), labels=[1,0,0,0], size=11,\

linewidth=0.1)
m.drawmeridians(np.arange(0, 360, 30.),labels=[0,1,0,1], size=11, \
linewidth=0.1)
m.drawcoastlines()
m.drawmapboundary()

Out[33]: <matplotlib.patches.Rectangle at 0x7f6ca1a66590>

0.1.5 Cmo eliminar los datos NaN de una matriz o un arreglo?

In [34]: print np.mean(Mapa)

No_NaN = Mapa[np.isfinite(Mapa)]
print np.mean(No_NaN)
print np.median(No_NaN)

4
nan
0.665697
0.43

0.1.6 Construccin de los vectores para la caracterizacin de la serie de datos

title

5
title

6
In [35]: import scipy.stats

NoNaN = []
Mean = []
Median= []
desv = []
perc_25 = []
perc_75 = []
kurt = []

for i in range(len(time)):
Mapa = tempanomaly[i,:,:]
Mapa_NoNaN = Mapa[np.isfinite(Mapa)]

NoNaN.append(len(Mapa_NoNaN))

Mean.append(np.mean(Mapa_NoNaN))
Median.append(np.median(Mapa_NoNaN))
desv.append(np.std(Mapa_NoNaN))
perc_25.append(np.percentile(Mapa_NoNaN,25))
perc_75.append(np.percentile(Mapa_NoNaN,75))
kurt.append(scipy.stats.kurtosis(Mapa_NoNaN))

NoNaN = np.array(NoNaN)
Mean = np.array(Mean)
Median = np.array(Median)
desv = np.array(desv)
perc_25 = np.array(perc_25)

7
 perc_75 = np.array(perc_75)
kurt = np.array(kurt)

yk = (perc_25 - 2*(Median) + perc_75) / (perc_75- perc_25)

0.1.7 Graficacin
In [38]: plt.figure(figsize=[6,7])
plt.plot(fechas,Mean, 'g', label = 'Media')
plt.plot(fechas,Median, 'b', label = 'Mediana')
#plt.plot(fechas,Mean + desv,'r',label=u'Desviacin Estandar')
#plt.plot(fechas,Mean - desv,'r')

plt.plot(fechas, perc_25, color='r',lw=0.4)

plt.plot(fechas, perc_75, color='r',lw=0.4)
plt.fill_between(fechas, Mean - desv,Mean + desv, alpha=0.2)
plt.legend(loc = 2)
plt.xlabel('Tiempo', fontsize=13)
plt.ylabel(ur'Anomalas Temp $\circ$ C', fontsize=13)
plt.title(u'Anomalas globales de Temperatura', fontsize=13)

Out[38]: <matplotlib.text.Text at 0x7f6c9b960950>

8
0.1.8 Graficacin de la serie de temperatura de un pixel
In [21]: Serie=tempanomaly[:,45,50]

plt.close('all')
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(fechas,Serie,'m',lw=1, label="Anomalias Temperatura")
plt.title('Serie de Anomalias Temperatura (GISTEMP)')
ax.set_xlabel('Tiempo')

9
 ax.set_ylabel('Anomalias Temperatura')

plt.show()

1 PROBABILIDAD
Como se mencion brevemente la primea reunin, en anlisis de datos ambientales trataremos de
hacer uso de herramientas estadsticas descriptivas y de inferencia con el objetivo de encontrar o
refutar la existencia determinada evidencia de inters en datos experimentales, comprobar hipte-
sis o, en general, mejorar el entendimiendo de los procesos fsicos que dieron lugar a determinados
registros experimentales. Un buen comienzo para el curso de anlisis es el recordar y hacer de los
conceptos fundamentales de Probabilidad (Breve historia del comienzo de la Teora de la Proba-
biliad: Siglo XVII, Pascal,Fermat, correspondencia sobre Juegos de Azar). Una de las preguntas
ms recurrentes de la sociedad cuando se habla de fenmenos natuarles est relacionada con la
probabilidad de ocurrencia de un determinado fenmeno. Desde el punto de vista de la gestin de
recursos, gestin de riesgos, desarrollo sostenible, desarrollo econmico y social de determinada
regin, es indispensable tener estimaciones sobre la probabilidad de ocurrencia, por ejemplo, de
un sismo de magnitud considerable, o dada su relevancia en Colombia, cual es la probabilidad
de desarrollo de un evento El Nio, cual es la probabilidad de ocurrencia de un evento de precip-
itacin intenso que afecte a una comunidad. La mirada clsica del concepto de probabilidad, y tal
vez la que permite mayor interpretacin en el campo de las aplicaciones incluyendo el anlisis de
datos en geociencias, est relacionado con la frecuencia de ocurrencia de un evento determinado.

10
En principio, y de manera resumida, se puede definir la probabilidad como la frecuencia relativa
de largo plazo y se puede escribir de la siguiente manera

2 n a  o
P r  P r {E}  0
 
n
cuando n
(Discusin: La teora de la probabilidad puede ser pensada a la vez como una herramienta de-
scriptiva como una herramientade inferencia.. Teora de Bayes)

title

3 PROBABILIDAD Y CONJUNTOS
Espacio de muestreo: S
Eventos: No precipitacin, precipitacin lquida, precipitacin congelada. . . (Cubren todas las
posibilidades: Evetos mutuamente excluyentes, colectivamente exhaustivos MECE)
Axiomas:
Probabilidad de todos los posibles eventos es 1.
Probabilidad de varios eventos mutuamente exlcuyentes es la suma de las probabilidades de
ocurrencia individual

11
title

title

12
 title

4
P r {E2 } < P r {E1 }

5
P r {E1 } < P r {E2 }

6
P r {E2 U E1 }??

7
P r {E2 } + P r {E1 }

8
P r {E2 U E1 }??

9
P r {E2 } + P r {E1 } P r {E2 E1 }

13
title

title

14
title

title

15
 title

9.0.1 Histogramas - Funciones de densidad de probabilidad

In [45]: # Seleccin de un mapa
a = np.where(fechas == datetime.datetime(1992,12,15))
print a
print a[0][0]

Mapa = tempanomaly[a[0][0],:,:]

# Se eliminan los NaN

SerieDeMapa = Mapa[np.isfinite(Mapa)]

print SerieDeMapa
print SerieDeMapa.shape

# Creacin del Histograma (Se pueden usar dos funciones, una del paquete
#numpy y otra de matplotlib. Se recomienda usar np.histogram)

num_bins=30
n1, bins1, patches1 = plt.hist(SerieDeMapa, num_bins, normed=True, \
facecolor='red', alpha=0.5)

ValueProb= 0.

fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(bins1[0:len(bins1)-1],n1/sum(n1),'r',lw=2, markersize=3, \

16
 label='Anomalias T') #f1

temp= n1/sum(n1)
plt.title('Pr{Anomalia Temperatura>='+np.str(ValueProb)+'}='+ \
np.str(len(SerieDeMapa[SerieDeMapa>=ValueProb])\
/np.float(len(SerieDeMapa))))

donde=[bins1[0:len(bins1)-1] >= ValueProb]

plt.fill_between((bins1[0:len(bins1)-1])[donde],(n1/sum(n1))[donde], \
0.0001, alpha=0.5,facecolor='gray', edgecolor='gray')

plt.legend(prop={'size':10})

plt.savefig('Histograma02.pdf')

(array([1355]),)
1355
[-1.8499999 -1.8499999 -1.8499999 ..., -0.81 -0.81 -0.81 ]
(16052,)

17
In [46]: num_bins = 20
hist,bins= np.histogram(SerieDeMapa, bins=15)
hist =hist.astype(float)
pdf= hist/np.sum(hist)

print len(bins)
print len(pdf)
print pdf
print bins

16
15
[ 0.00274109 0.00865936 0.02093197 0.08011463 0.09973835 0.1716297
0.30837279 0.18477448 0.05320209 0.02155495 0.01937453 0.01065288
0.00778719 0.00672813 0.00373785]
[-3.66999984 -3.07799985 -2.48599985 -1.89399986 -1.30199987 -0.70999988
-0.11799989 0.4740001 1.0660001 1.65800009 2.25000008 2.84200007
3.43400006 4.02600005 4.61800005 5.21000004]

In [59]: fig = plt.figure(figsize=[9,6])

ax = fig.add_subplot(111)
temp=pdf/np.sum(pdf)
print temp

18
bincenters=(bins[1:]+bins[:-1])/2
plt.plot(bincenters,pdf,'o-',color='b',lw=2, markersize=10)

from scipy.interpolate import interp1d

fdi= interp1d(bincenters,pdf)

plt.plot([np.median(SerieDeMapa),np.median(SerieDeMapa)] \
,[0,fdi(np.median(SerieDeMapa))],color='g',lw=3)

donde=np.where((bincenters>=np.percentile(SerieDeMapa,10)) & \
(bincenters<=np.percentile(SerieDeMapa,50)))

print np.concatenate((np.array([np.percentile(SerieDeMapa,10)]),\
bincenters[donde],\
(np.array([np.percentile(SerieDeMapa,50)]))))

xfill=np.concatenate((np.array([np.percentile(SerieDeMapa,10)]),\
bincenters[donde],\
(np.array([np.percentile(SerieDeMapa,50)]))))

yfill=np.concatenate((np.array([fdi(np.percentile(SerieDeMapa,10))]),\
pdf[donde],\
(np.array([fdi(np.percentile(SerieDeMapa,50))]))))

plt.fill_between(xfill,yfill, 0.0001, alpha=0.5,facecolor='gray', \

edgecolor='gray')

donde2=np.where((bincenters>=np.percentile(SerieDeMapa,50)) & \
(bincenters<=np.percentile(SerieDeMapa,90)))

xfill=np.concatenate((np.array([np.percentile(SerieDeMapa,50)]),\
bincenters[donde2],\
(np.array([np.percentile(SerieDeMapa,90)]))))

yfill=np.concatenate((np.array([fdi(np.percentile(SerieDeMapa,50))]),\
pdf[donde2],\
(np.array([fdi(np.percentile(SerieDeMapa,90))]))))

plt.fill_between(xfill,yfill, \
0.0001, alpha=0.5,\
facecolor='m', edgecolor='m')

19
 plt.title('Histograma Serie Mapa')

plt.savefig('Histograma03.pdf')

[ 0.00274109 0.00865936 0.02093197 0.08011463 0.09973835 0.1716297

0.30837279 0.18477448 0.05320209 0.02155495 0.01937453 0.01065288
0.00778719 0.00672813 0.00373785]
[-1.49000001 -1.00599988 -0.41399988 0.11 ]

In [36]: np.percentile(SerieDeMapa, 25)

Out[36]: -0.4699999988079071

10 PROBABILIDAD CONDICIONAL
Probabilidad de un evento dado que otro ocurri u ocurrir
P r {E2 |E1 } = Probabilidad de E2 dado que E1 ocurri u ocurrir.

20
11
P r {E2 E1 }
P r {E2 |E1 } =
P r {E1 }

12 CUANTILES

13
qp

Valor que excede la fraccin de datos correspondiente a la probabilidad p

Cmo programar esto?? Cual es la manera mas fcil?? Como lo haria manualmente?
Datos ordenados, datos ranqueados

14 
{x1 , x2 , x3 , ...., xn } x(1) , x(2) , x(3) , ...., x(n)

21