Dans les articles précédents, j'ai abordé certaines manières de traiter Données de l'Agence météorologique, mais il existe une autre méthode d'acquisition de données (état actuel). , Cette méthode est la plus efficace), et je vais alléger l'analyse des données acquises.

Articles passés

・ Météorologie x Python ~ Acquisition automatique des données de localisation AMeDAS ~ https://qiita.com/OSAKO/items/264c77b70843045bc12b ・ Météorologie x Python ~ Acquisition automatique des données de points AMeDAS (édition supplémentaire) ~ https://qiita.com/OSAKO/items/505ecee67df424963e53 ・ Météorologie x Ruby ~ Grattage de rubis à l'aide de Mechanize ~ https://qiita.com/OSAKO/items/3c1cac0b5448be9ab243

1. Ramper

▶ Je veux obtenir les valeurs numériques ou les chaînes de caractères incorporées dans le format tabulaire suivant à la fois. (Exemple: données du 1er janvier 2019 à Tokyo) https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=44&block_no=47662&year=2019&month=01&day=1&view=

▶ Tout d'abord, implémentez un robot d'exploration pour obtenir toutes les URL de plus de 1000 points et les lister. https://www.data.jma.go.jp/obd/stats/etrn/index.php?prec_no=44&block_no=47662&year=&month=&day=&view= ↑ Exemple d'url que vous souhaitez obtenir ( ** prec_no et block_no les identifiants diffèrent selon l'emplacement ** </ font>)

procédure

① Tout d'abord, depuis [cette url](https://www.data.jma.go.jp/obd/stats/etrn/select/prefecture00.php?prec_no=0block_no Tout le monde & year = & month = & day = & view =) Je veux une URL pour accéder, donc je recherche l'attribut href sous la balise area et je le liste. (Méthode Get_area_link dans le code ci-dessous) (2) En utilisant l'url de chaque zone listée ci-dessus comme clé, si vous recherchez l'attribut href de chaque point de la même manière, vous pouvez obtenir l'url cible avec id dans prec_no et block_no. (Méthode Get_station_link dans le code ci-dessous) ③ Après cela, formatez un peu l'url répertoriée et enregistrez-la au format csv. (Méthode Data_arange dans le code ci-dessous)

programme

`get_amedas_station_url.py`


# -*- coding: utf-8 -*-
import pandas as pd
import urllib.request
from bs4 import BeautifulSoup

class Get_amedas_station:
  def __init__(self):
    url = 'https://www.data.jma.go.jp/obd/stats/etrn/select/prefecture00.php?prec_no=&block_no=&year=&month=&day=&view='
    html = urllib.request.urlopen(url)
    self.soup = BeautifulSoup(html, 'html.parser')
  
  def get_area_link(self):
    elements = self.soup.find_all('area')
    self.area_list = [element['alt'] for element in elements]
    self.area_link_list = [element['href'] for element in elements]

  def get_station_link(self):
    out = pd.DataFrame(columns=['station','url','area'])
    for area, area_link in zip(self.area_list, self.area_link_list):
      url = 'https://www.data.jma.go.jp/obd/stats/etrn/select/'+ area_link
      html = urllib.request.urlopen(url)
      soup = BeautifulSoup(html, 'html.parser')
      elements = soup.find_all('area')
      station_list = [element['alt'] for element in elements]
      station_link_list = [element['href'].strip('../') for element in elements]
      df1 = pd.DataFrame(station_list,columns=['station'])
      df2 = pd.DataFrame(station_link_list,columns=['url'])
      df = pd.concat([df1, df2],axis=1).assign(area=area)
      out = pd.concat([out,df])
      print(area)
    self.out = out
  
  def data_arange(self):
    out = self.out[~self.out.duplicated()].assign(append='https://www.data.jma.go.jp/obd/stats/etrn/')
    out['amedas_url'] = out['append'] + out['url']
    out = out.loc[:,['area','station','amedas_url']]
    out.to_csv('amedas_url_list.csv',index=None, encoding='SJIS')

if __name__ == '__main__':
  amedas = Get_amedas_station()
  amedas.get_area_link()
  amedas.get_station_link()
  amedas.data_arange()

2. Grattage

▶ En ce qui concerne l'acquisition de données AMeDAS, qui est le sujet principal, cette fois j'ai écrit un programme qui prend en charge différentes échelles de temps de données de 1 jour, 1 heure et 10 minutes (cependant, la quantité de code est relativement importante et cela peut être un gaspillage, donc c'est nécessaire. Amélioration) ▶ En guise de mise en garde, l'élément d'observation diffère selon le point AMeDAS sélectionné. ( ** Le point sélectionné est différent selon qu'il s'agit d'une station météorologique ou non ** </ font>)

Par exemple, comparons les événements observés entre Tokyo et Hachioji sur chaque échelle de temps.

** ・ Données 1 jour ** Tokyo Hachioji ** ・ 1 heure de données ** Tokyo Hachioji ** ・ 10 minutes de données ** Tokyo Hachioji

programme

`get_amedas_data.py`


# -*- coding: utf-8 -*-
import re
import time
import numpy as np
import pandas as pd
import urllib.request
from bs4 import BeautifulSoup
from datetime import timedelta
from datetime import datetime as dt
from dateutil.relativedelta import relativedelta

#Une méthode pour rechercher et lister les données à acquérir à partir de l'url de chaque point AMeDAS
def search_data(url):
  html = urllib.request.urlopen(url)
  time.sleep(1)
  soup = BeautifulSoup(html, 'html.parser')
  element = soup.find_all('tr', attrs={'class':'mtx', 'style':'text-align:right;'})
  out = [list(map(lambda x: x.text, ele)) for ele in element]
  return out

class Get_amedas_data:
  def __init__(self,area_name,station_name):
    #Faites référence à l'url cible de la zone et du point spécifiés
    self.st_name = station_name
    self.a_name = area_name
    amedas_url_list = pd.read_csv('amedas_url_list.csv',encoding='SJIS')
    df = amedas_url_list[(amedas_url_list['area']==self.a_name) & (amedas_url_list['station']==self.st_name)]
    amedas_url = df.iat[0,2]
    #Obtenir des nombres avec des expressions régulières à partir de l'URL (prec_non et bloquer_aucun identifiant requis)
    pattern=r'([+-]?[0-9]+\.?[0-9]*)'
    id_list=re.findall(pattern, amedas_url)
    self.pre_id = id_list[0]
    self.s_id = id_list[1]

  #Étant donné que les données d'un jour peuvent être acquises collectivement en tant que données de janvier, spécifiez le mois de début et le mois de fin.
  def set_date1(self, startmonth, endmonth):
    self.start = startmonth
    self.end = endmonth
    strdt = dt.strptime(self.start, '%Y%m')
    enddt = dt.strptime(self.end, '%Y%m')
    months_num = (enddt.year - strdt.year)*12 + enddt.month - strdt.month + 1
    #Lister les mois du mois de début au mois de fin
    self.datelist = map(lambda x, y=strdt: y + relativedelta(months=x), range(months_num))
  
  #Pour les données de 1 heure ou 10 minutes, spécifiez la date de début et la date de fin
  def set_date2(self,startdate,enddate):
    self.start = startdate
    self.end = enddate
    strdt = dt.strptime(self.start, '%Y%m%d')
    enddt = dt.strptime(self.end, '%Y%m%d')
    days_num = (enddt - strdt).days + 1
    #Liste des dates de la date de début à la date de fin
    self.datelist = map(lambda x, y=strdt: y + timedelta(days=x), range(days_num))

  #Créez à l'avance un bloc de données vide. Il y a plus d'éléments qui peuvent être acquis au point où se trouve l'observatoire météorologique.
  #Spécifiez l'échelle de temps que vous souhaitez obtenir (type)
  def dl_data(self, type):
    #Trame de données vide pour les données quotidiennes
    data1  = pd.DataFrame(columns=['Année mois','journée','Pression locale moyenne','Pression moyenne au niveau de la mer','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Humidité moyenne','Humidité minimale','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','journée照時間','chute de neige','Neige la plus profonde','Aperçu météo (midi)','Aperçu météo (nuit)'])
    data1_ = pd.DataFrame(columns=['Année mois','journée','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','La plupart des directions du vent','journée照時間','chute de neige','Neige la plus profonde'])
    #Trame de données vide pour les données horaires
    data2  = pd.DataFrame(columns=['Date','Temps','Pression locale','Pression au niveau de la mer','Précipitation','Température','Température du point de rosée','La pression de la vapeur','Humidité','vitesse du vent','Direction du vent','日照Temps間','Rayonnement solaire total','chute de neige','La couverture de neige','Météo','Volume du cloud','Visibilité'])
    data2_ = pd.DataFrame(columns=['Date','Temps','Précipitation','Température','vitesse du vent','Direction du vent','日照Temps間','chute de neige','La couverture de neige'])
    #Trame de données vide pour 10 min de données
    data3  = pd.DataFrame(columns=['Date','Heure et minute','Pression locale','Pression au niveau de la mer','Précipitation','Température','Humidité relative','Vitesse moyenne du vent','Direction moyenne du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','Heure du soleil'])
    data3_ = pd.DataFrame(columns=['Date','Heure et minute','Précipitation','Température','Vitesse moyenne du vent','Direction moyenne du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','Heure du soleil'])

    #Créez un bloc de données en le rejoignant verticalement lors de l'acquisition de données tout en faisant pivoter le mois ou la liste de dates répertoriés
    for dt in self.datelist:
      d = dt.strftime("%Y%m%d")
      yyyy = d[0:4]
      mm = d[4:6]
      dd = d[6:8]

      if type=='daily':
        #Bloquer au point où se trouve l'observatoire météorologique_non est un nombre à 5 chiffres
        if len(self.s_id) == 5:
          pattern = 's1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          df = (pd.DataFrame(out, columns=['journée','Pression locale moyenne','Pression moyenne au niveau de la mer','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Humidité moyenne','Humidité minimale','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','journée照時間','chute de neige','Neige la plus profonde','Aperçu météo (midi)','Aperçu météo (nuit)'])).assign(Année mois=f'{yyyy}{mm}')
          df = df.loc[:,['Année mois','journée','Pression locale moyenne','Pression moyenne au niveau de la mer','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Humidité moyenne','Humidité minimale','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','journée照時間','chute de neige','Neige la plus profonde','Aperçu météo (midi)','Aperçu météo (nuit)']]
          data1 = pd.concat([data1, df])
          data1.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')
        else:
          pattern = 'a1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          df = (pd.DataFrame(out, columns=['journée','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','La plupart des directions du vent','journée照時間','chute de neige','Neige la plus profonde'])).assign(Année mois=f'{yyyy}{mm}')
          df = df.loc[:,['Année mois','journée','journée降水量','Jusqu'à 1 heure de précipitations','Jusqu'à 10 minutes de précipitations','Température moyenne','Température la plus élevée','Température la plus basse','Vitesse moyenne du vent','Vitesse maximale du vent','Direction maximale du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','La plupart des directions du vent','journée照時間','chute de neige','Neige la plus profonde']]
          data1_ = pd.concat([data1_, df])
          data1_.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')

      elif type=='hourly':
        if len(self.s_id) == 5:
          pattern = 's1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          #Colonne de date ajoutée (24)
          date = pd.DataFrame((np.full([24,1], f'{yyyy}{mm}{dd}')),columns=['Date'])
          df = pd.DataFrame(out, columns=['Temps','Pression locale','Pression au niveau de la mer','Précipitation','Température','Température du point de rosée','La pression de la vapeur','Humidité','vitesse du vent','Direction du vent','日照Temps間','Rayonnement solaire total','chute de neige','La couverture de neige','Météo','Volume du cloud','Visibilité'])
          df = pd.concat([date, df],axis=1)
          data2 = pd.concat([data2, df])
          data2.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')
        else:
          pattern = 'a1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          date = pd.DataFrame((np.full([24,1], f'{yyyy}{mm}{dd}')),columns=['Date'])
          df = pd.DataFrame(out, columns=['Temps','Précipitation','Température','vitesse du vent','Direction du vent','日照Temps間','chute de neige','La couverture de neige'])
          df = pd.concat([date, df],axis=1)
          data2_ = pd.concat([data2_, df])
          data2_.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')

      elif type=='10min':
        if len(self.s_id) == 5:
          pattern = 's1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          #Colonne de date ajoutée (6 x 24)
          date = pd.DataFrame((np.full([144,1], f'{yyyy}{mm}{dd}')),columns=['Date'])
          df = pd.DataFrame(out, columns=['Heure et minute','Pression locale','Pression au niveau de la mer','Précipitation','Température','Humidité relative','Vitesse moyenne du vent','Direction moyenne du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','Heure du soleil'])
          df = pd.concat([date, df],axis=1)
          data3 = pd.concat([data3, df])
          data3.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')
        else:
          pattern = 'a1'
          url = f'https://www.data.jma.go.jp/obd/stats/etrn/view/{type}_{pattern}.php?prec_no={self.pre_id}&block_no={self.s_id}&year={yyyy}&month={mm}&day={dd}&view=p1'
          out = search_data(url)
          date = pd.DataFrame((np.full([144,1], f'{yyyy}{mm}{dd}')),columns=['Date'])
          df = pd.DataFrame(out, columns=['Heure et minute','Précipitation','Température','Vitesse moyenne du vent','Direction moyenne du vent','Vitesse du vent instantanée maximale','Direction instantanée maximale du vent','Heure du soleil'])
          df = pd.concat([date, df],axis=1)
          data3_ = pd.concat([data3_, df])
          data3_.to_csv(f'{self.a_name}_{self.st_name}_{self.start}-{self.end}_{type}.csv',index=None, encoding='SJIS')

      print(f'{self.a_name}_{self.st_name}_{yyyy}-{mm}-{dd}_{type}')
    print(f'{self.a_name}_{self.st_name}Dans{self.start}〜{self.end}de{type}J'ai téléchargé les données.')


if __name__ == '__main__':
  a_name = input('Entrez la région que vous souhaitez télécharger:')
  st_name = input('Entrez l'emplacement que vous souhaitez télécharger:')
  amedas = Get_amedas_data(a_name,st_name)
  type = input('Sélectionnez une échelle de temps (quotidienne ou horaire ou 10 min):')
  if type == 'daily':
    start= input('Mois de début que vous souhaitez obtenir(yyyymm)Entrez s'il vous plait:')
    end = input('Fin du mois que vous souhaitez obtenir(yyyymm)Entrez s'il vous plait:')
    amedas.set_date1(start,end)
  else:
    start= input('Date de début que vous souhaitez obtenir(Veuillez saisir aaaammjj:')
    end = input('Date de fin que vous souhaitez obtenir(yyyymmdd)Entrez s'il vous plait:')
    amedas.set_date2(start,end)
  
  amedas.dl_data(type)

Région: Tokyo Point: Tokyo Échelle de temps: 10min Date de début: 20190401 Date de fin: 20190731

Je vais essayer de télécharger les données sous la condition.

Le contenu du fichier "Tokyo_Tokyo_20190401-20190731_10min.csv" enregistré ressemble à ceci ![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/354874/10817920-dcbe-ba75-a2f8-34071b89d19d.png) # 3. Analyse du spectre ▶ Comme j'ai acquis des données détaillées avec une résolution temporelle de 10 minutes, j'analyserai légèrement les données de la série chronologique. ▶ Cette fois, le spectre de puissance a été calculé en appliquant une analyse spectrale à la fluctuation de température pendant la période acquise. Pour ceux qui veulent comprendre les sentiments de la transformation de Fourier, [cette vidéo](https://www.youtube.com/watch?v=bjBZEKdlLD0&t=339s) peut être facile à comprendre.

`fft_temp.py`


# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from numpy.fft import fftn, ifftn, fftfreq

class FFT:
  def __init__(self):
    df = pd.read_csv('Tokyo_Tokyo_20190401-20190731_10min.csv', encoding='SJIS').replace(['--','×'],-99)
    df['Température'] = df['Température'].astype(float)
    y = df['Température'].replace(-99, (df['Température'] > -99).mean()).values
    y = (y - np.mean(y)) * np.hamming(len(y))
    '''
Fréquence d'échantillonnage fs → Nombre de données échantillonnées par seconde d dans le deuxième argument=1.0/fs
Étant donné que la fréquence d'échantillonnage est cette fois des données d'intervalle de 10 minutes, 1/600
    '''
    self.z = np.fft.fft(y)
    self.freq = fftfreq(len(y), d=1.0/(1/600))


  def plot(self, n_samples):
    fig, axes = plt.subplots(figsize=(8, 4), ncols=2, sharey=True)
    ax = axes[0]
    ax.plot(self.freq[1:int(n_samples/2)], abs(self.z[1:int(n_samples/2)]))
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_ylabel('Power')
    ax.set_xlabel('Fréquence (Hz)')
    ax = axes[1]
    ax.plot((1 / self.freq[1:int(n_samples / 2)])/(60*60), abs(self.z[1:int(n_samples / 2)]))
    ax.set_yscale('log')
    ax.set_xscale('log')
    ax.set_xlabel('Temps d'un cycle)')
    plt.show()
  
  def periodic_characteristics(self, n_samples):
    fft_pow_df = pd.DataFrame([(1 / self.freq[1:int(n_samples / 2)])/(60*60), np.log10(abs(self.z[1:int(n_samples / 2)]))], index=['Temps d'un cycle)', 'log10_power']).T
    fft_pow_df = fft_pow_df.sort_values('log10_power', ascending=False).head(10).reset_index()
    print(fft_pow_df.loc[:, ['Temps d'un cycle)', 'log10_power']])

if __name__ == '__main__':
  fft = FFT()
  fft.periodic_characteristics(512)
  fft.plot(512)

Le cycle de 122 jours, qui a le spectre de puissance le plus large, n'est pas tout à fait correct, mais comme prévu, le cycle de 23,24 heures était également exceptionnel, alors je me demande s'il convient. (Parce qu'il a les caractéristiques physiques d'un cycle quotidien que la température monte du matin au midi et redescend de nuit) Je vous serais reconnaissant si vous pouviez signaler des erreurs.

c'est tout! !! !!

Météorologie x Python ~ De l'acquisition de données météorologiques à l'analyse spectrale ~

1. Ramper

procédure

programme

get_amedas_station_url.py

2. Grattage

programme

get_amedas_data.py

fft_temp.py

`get_amedas_station_url.py`

`get_amedas_data.py`

`fft_temp.py`