TALLER 7 - Python

Geoinformática
Pedro Joaquin Perilla Vargas
Métodos de Fourier

February 7, 2020

1. Lea los datos de sunspot_year.dat, en la segunda columna. La primera columna representa

el año de medición. Grafique. Calcule el espectro utilizando Periodograma con y sin
Hanning taper amp_spec y amp_spec_taper. Utilice como dt = 1 año. Grafique Periodo vs
Espectro, entre periodos de 1 y 20 años. Note que periodo es 1/freq.

Solución:

Lectura de los datos:

[1]: import numpy as np

sunspotdata='sunspot_year.dat'

# Lee el archivo
data = np.loadtxt(sunspotdata)
year = data[:,0]
ms = data[:,1]
print(year)
print(ms)

[1700. 1701. 1702. 1703. 1704. 1705. 1706. 1707. 1708. 1709. 1710. 1711.
1712. 1713. 1714. 1715. 1716. 1717. 1718. 1719. 1720. 1721. 1722. 1723.
1724. 1725. 1726. 1727. 1728. 1729. 1730. 1731. 1732. 1733. 1734. 1735.
1736. 1737. 1738. 1739. 1740. 1741. 1742. 1743. 1744. 1745. 1746. 1747.
1748. 1749. 1750. 1751. 1752. 1753. 1754. 1755. 1756. 1757. 1758. 1759.
1760. 1761. 1762. 1763. 1764. 1765. 1766. 1767. 1768. 1769. 1770. 1771.
1772. 1773. 1774. 1775. 1776. 1777. 1778. 1779. 1780. 1781. 1782. 1783.
1784. 1785. 1786. 1787. 1788. 1789. 1790. 1791. 1792. 1793. 1794. 1795.
1796. 1797. 1798. 1799. 1800. 1801. 1802. 1803. 1804. 1805. 1806. 1807.
1808. 1809. 1810. 1811. 1812. 1813. 1814. 1815. 1816. 1817. 1818. 1819.
1820. 1821. 1822. 1823. 1824. 1825. 1826. 1827. 1828. 1829. 1830. 1831.
1832. 1833. 1834. 1835. 1836. 1837. 1838. 1839. 1840. 1841. 1842. 1843.
1844. 1845. 1846. 1847. 1848. 1849. 1850. 1851. 1852. 1853. 1854. 1855.
1856. 1857. 1858. 1859. 1860. 1861. 1862. 1863. 1864. 1865. 1866. 1867.

1
 1868. 1869. 1870. 1871. 1872. 1873. 1874. 1875. 1876. 1877. 1878. 1879.
1880. 1881. 1882. 1883. 1884. 1885. 1886. 1887. 1888. 1889. 1890. 1891.
1892. 1893. 1894. 1895. 1896. 1897. 1898. 1899. 1900. 1901. 1902. 1903.
1904. 1905. 1906. 1907. 1908. 1909. 1910. 1911. 1912. 1913. 1914. 1915.
1916. 1917. 1918. 1919. 1920. 1921. 1922. 1923. 1924. 1925. 1926. 1927.
1928. 1929. 1930. 1931. 1932. 1933. 1934. 1935. 1936. 1937. 1938. 1939.
1940. 1941. 1942. 1943. 1944. 1945. 1946. 1947. 1948. 1949. 1950. 1951.
1952. 1953. 1954. 1955. 1956. 1957. 1958. 1959. 1960. 1961. 1962. 1963.
1964. 1965. 1966. 1967. 1968. 1969. 1970. 1971. 1972. 1973. 1974. 1975.
1976. 1977. 1978. 1979. 1980. 1981. 1982. 1983. 1984. 1985. 1986. 1987.
1988. 1989. 1990. 1991. 1992. 1993. 1994. 1995. 1996. 1997. 1998. 1999.
2000. 2001. 2002. 2003. 2004. 2005. 2006. 2007.]
[ 5. 11. 16. 23. 36. 58. 29. 20. 10. 8. 3. 0.
0. 2. 11. 27. 47. 63. 60. 39. 28. 26. 22. 11.
21. 40. 78. 122. 103. 73. 47. 35. 11. 5. 16. 34.
70. 81. 111. 101. 73. 40. 20. 16. 5. 11. 22. 40.
60. 80.9 83.4 47.7 47.8 30.7 12.2 9.6 10.2 32.4 47.6 54.
62.9 85.9 61.2 45.1 36.4 20.9 11.4 37.8 69.8 106.1 100.8 81.6
66.5 34.8 30.6 7. 19.8 92.5 154.4 125.9 84.8 68.1 38.5 22.8
10.2 24.1 82.9 132. 130.9 118.1 89.9 66.6 60. 46.9 41. 21.3
16. 6.4 4.1 6.8 14.5 34. 45. 43.1 47.5 42.2 28.1 10.1
8.1 2.5 0. 1.4 5. 12.2 13.9 35.4 45.8 41.1 30.1 23.9
15.6 6.6 4. 1.8 8.5 16.6 36.3 49.6 64.2 67. 70.9 47.8
27.5 8.5 13.2 56.9 121.5 138.3 103.2 85.7 64.6 36.7 24.2 10.7
15. 40.1 61.5 98.5 124.7 96.3 66.6 64.5 54.1 39. 20.6 6.7
4.3 22.7 54.8 93.8 95.8 77.2 59.1 44. 47. 30.5 16.3 7.3
37.6 74. 139. 111.2 101.6 66.2 44.7 17. 11.3 12.4 3.4 6.
32.3 54.3 59.7 63.7 63.5 52.2 25.4 13.1 6.8 6.3 7.1 35.6
73. 85.1 78. 64. 41.8 26.2 26.7 12.1 9.5 2.7 5. 24.4
42. 63.5 53.8 62. 48.5 43.9 18.6 5.7 3.6 1.4 9.6 47.4
57.1 103.9 80.6 63.6 37.6 26.1 14.2 5.8 16.7 44.3 63.9 69.
77.8 64.9 35.7 21.2 11.1 5.7 8.7 36.1 79.7 114.4 109.6 88.8
67.8 47.5 30.6 16.3 9.6 33.2 92.6 151.6 136.3 134.7 83.9 69.4
31.5 13.9 4.4 38. 141.7 190.2 184.8 159. 112.3 53.9 37.6 27.9
10.2 15.1 47. 93.8 105.9 105.5 104.5 66.6 68.9 38. 34.5 15.5
12.6 27.5 92.5 155.4 154.6 140.4 115.9 66.6 45.9 17.9 13.4 29.4
100.2 157.6 142.6 145.7 94.3 54.6 29.9 17.5 8.6 21.5 64.3 93.3
119.6 111. 104. 63.7 40.4 29.8 15.2 7.5]

Gráfica de los datos de sunspot_year.dat

[2]: import matplotlib

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(year,ms, color='red')
ax.set(xlabel='Year', ylabel='Sunspot', title='Sunspot per year')
ax.grid()

2
 plt.show()

[3]: import amp_spec

import math
import numpy as np
import matplotlib.pyplot as plt

# Load file
fname = 'sunspot_year.dat'
data = np.loadtxt(fname,usecols=[1])
npts = np.size(data)

# Serie de tiempo
dt = 1 # año
t = np.arange(npts)*dt

# Calcular el espectro de amplitudes

freq, spec = amp_spec.amp_spec(data,dt)
Per = np.zeros(len(freq))
Per[1:] = 1./freq[1:]
#
# Calcular el espectro de amplitudes Hanning
freq1, spec1 = amp_spec.amp_spec_taper(data,dt)

3
 Per1 = np.zeros(len(freq))
Per1[1:] = 1./freq[1:]

# Figure de los datos

fig = plt.figure()
ax1 = fig.add_subplot(2,1,1)
ax1.plot(Per[1:],spec[1:],'r')
ax1.plot(Per1[1:],spec1[1:])
ax1.set_xlim(0, 20)
ax1.set_xlabel('Periodo (años)')
plt.show()

Conclusión: En el análisis de los datos en un periodo entre 1 y 20 años nos damos cuenta de
que sin el filtro de hanning se obtienen picos altos que representan información que no es
real, al aplicar el filtro se suaviza la señal y se ve de manera más adecuada la información.

2. En meteorología y climatología se estudia el régimen de temperaturas de lluvias y precip-

itación. En particular en la región andina, el régimen bimodal es común, pero pueden existir
lugares con regímenes diferentes (unimodal por ejemplo)>Esto se observa en los histogra-
mas de precipitación o temperatura.
Hacer histogramas de varias estaciones para mostrar estos regímenes. Verificar el compor-
tamiento cíclico con los espectros de amplitudes para esas estaciones. Interprete.

Solución:

Gráficas de los espectros de amplitudes amp_spec y amp_spec_taper para la estación

de interés, en este caso la 212.

[5]: import numpy as np

import matplotlib.pyplot as plt; plt.rcdefaults()
import cartopy.crs as ccrs

4
###---Inicio-Lectura de los archivos---###
Tfile = 'MON_T_CRU_19012015.csv'
Pfile = 'MON_P_CRU_19012015.csv'
sta_file = 'BASIN_CHARACTERISTICS.csv'

# Carga los archivos

basin = np.loadtxt(sta_file,skiprows=1,delimiter=',')
Temp = np.loadtxt(Tfile,skiprows=1,delimiter=',')
Prec = np.loadtxt(Pfile,skiprows=1,delimiter=',')

# Organiza
T = Temp[:,1:] # Filas, y 1ra columna
P = Prec[:,1:]

sta = basin[:,0] # Toma los datos de cada columna

lon = basin[:,1]
lat = basin[:,2]
area = basin[:,3]
elev = basin[:,4]
###---Fin-Lectura de los archivos---###

# Selección de estación
i = 212

#Crea el dt
dt=1/12 #Años

print(i,sta[i],lat[i],lon[i])

freq,Tspec =amp_spec.amp_spec(T[i,:],dt)
freq,Pspec =amp_spec.amp_spec(P[i,:],dt)
freq,Tspec1=amp_spec.amp_spec_taper(T[i,:],dt)
freq,Pspec1 =amp_spec.amp_spec_taper(P[i,:],dt)
Per = np.zeros(len(freq))
Per[i:]=1/freq[i:]

# Figure de los datos para Temperatura

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(Per[i:],Tspec[i:],'r')
ax1.plot(Per[i:],Tspec1[i:],'orange')
ax1.set_xlabel('Per')
ax1.set_ylabel('Tspec')

5
 ax2 = fig.add_subplot(212)
ax2.plot(Per[i:],Pspec[i:],'blue')
ax2.plot(Per[i:],Pspec1[i:],'black')
ax2.set_xlabel('Per')
ax2.set_ylabel('Pspec')

plt.show()

212 2907042.0 51.35 110.47

Creación de los histogramas correspondientes, en este caso para la estación 212.

[6]: import matplotlib.pyplot as plt
import numpy as np
month=np.arange(12)+1
i = 212
temp=T[i,:]
prec=P[i,:]

# Creación del mapa con la ubicación de las estacion o cuenca a nivel

,→global.

fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree())

6
ax.set_extent([-180, 180, -90, 90])
ax.coastlines()
plt.plot(lon,lat,'k^',markersize=10) # Convención para todas las estaciones
plt.plot(lon[i],lat[i],'r^',markersize=10) # Resaltando la estación de
,→interés

plt.show()

#Calculos histogramas
t_mean=np.zeros(12)
p_mean=np.zeros(12)
for imonth in month:
k=imonth-1
t_mean[k]=np.mean(temp[k::12])
p_mean[k]=np.mean(prec[k::12])

#Se crean las figuras de los histogramas

#Para los datos de temperatura

plt.subplot(211)
plt.bar(k, t_mean[k], color='red', edgecolor='black')
plt.ylabel('Temperatura')

#Para los datos de precipitación

plt.subplot(212)
plt.bar(k, p_mean[k], color='blue', edgecolor='black')
plt.ylabel('Precipitación')

7
 Gráficas para los espectros de amplitudes amp_spec y amp_spec_taper de la estacion
55
[7]: import numpy as np
import matplotlib.pyplot as plt; plt.rcdefaults()
import cartopy.crs as ccrs

###---Inicio-Lectura de los archivos---###

Tfile = 'MON_T_CRU_19012015.csv'
Pfile = 'MON_P_CRU_19012015.csv'
sta_file = 'BASIN_CHARACTERISTICS.csv'

# Carga los archivos

basin = np.loadtxt(sta_file,skiprows=1,delimiter=',')
Temp = np.loadtxt(Tfile,skiprows=1,delimiter=',')
Prec = np.loadtxt(Pfile,skiprows=1,delimiter=',')

# Organiza
T = Temp[:,1:] # Filas, y 1ra columna
P = Prec[:,1:]

sta = basin[:,0] # Toma los datos de cada columna

lon = basin[:,1]

8
 lat = basin[:,2]
area = basin[:,3]
elev = basin[:,4]
###---Fin-Lectura de los archivos---###

# Selección de estación
i = 55

#Crea el dt
dt=1/12 #Años

print(i,sta[i],lat[i],lon[i])

freq,Tspec =amp_spec.amp_spec(T[i,:],dt)
freq,Pspec =amp_spec.amp_spec(P[i,:],dt)
freq,Tspec1=amp_spec.amp_spec_taper(T[i,:],dt)
freq,Pspec1 =amp_spec.amp_spec_taper(P[i,:],dt)
Per = np.zeros(len(freq))
Per[i:]=1/freq[i:]

# Figure de los datos para Temperatura

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(Per[i:],Tspec[i:],'r')
ax1.plot(Per[i:],Tspec1[i:],'orange')
ax1.set_xlabel('Per')
ax1.set_ylabel('Tspec')

ax2 = fig.add_subplot(212)
ax2.plot(Per[i:],Pspec[i:],'blue')
ax2.plot(Per[i:],Pspec1[i:],'black')
ax2.set_xlabel('Per')
ax2.set_ylabel('Pspec')

plt.show()

55 1339500.0 3.52 11.5

9
 Creación de los histogramas correspondientes, en este caso para la estación 55.
[8]: import matplotlib.pyplot as plt
import numpy as np
month=np.arange(12)+1
i = 55
temp=T[i,:]
prec=P[i,:]

# Creación del mapa con la ubicación de las estacion o cuenca a

,→nivel global.
fig = plt.figure(figsize=(20,10))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree())
ax.set_extent([-180, 180, -90, 90])
ax.coastlines()
plt.plot(lon,lat,'k^',markersize=10) # Convención para todas las
,→estaciones

plt.plot(lon[i],lat[i],'r^',markersize=10) # Resaltando la estación

,→de interés

plt.show()

#Calculos histogramas

10
t_mean=np.zeros(12)
p_mean=np.zeros(12)
for imonth in month:
k=imonth-1
t_mean[k]=np.mean(temp[k::12])
p_mean[k]=np.mean(prec[k::12])

#Se crean las figuras de los histogramas

#Para los datos de temperatura

plt.subplot(211)
plt.bar(k, t_mean[k], color='red', edgecolor='black')
plt.ylabel('Temperatura')

#Para los datos de precipitación

plt.subplot(212)
plt.bar(k, p_mean[k], color='blue', edgecolor='black')
plt.ylabel('Precipitación')

11
Interpretación: Se realizó el ejercicio de los histrogramás y los espectros de amplitudes para
las estaciones 212 y 55, observando los histogramas correspondientes se llega a la conclusión
de que para la estación 212 se tiene un régimen unimodal lo cual se puede evidenciar por
el comportamiento de las barras del histograma de precipitación y para la estación 55 se
tiene un régimen bimodal igualmente deducido por el histograma de precipitación. Esta
información se contrastó con los espectros de amplitudes correspondiente verificándose el
comportamiento cíclico de las estaciones y además obteniendo una visión más detallada
para el periodo estudiado con los dos filtros utilizados.

12