Introduction du code de dessin pour les figures avec un certain degré de perfection des données météorologiques

Cet article est le 11ème jour du Calendrier de l'Avent étudiant LT 2019 - Qiita.

Bonjour. Je suis un étudiant diplômé spécialisé en sciences atmosphériques.

Cette fois, comme vous pouvez le voir dans l'article précédent, je me suis concentré sur le dessin de données météorologiques à l'aide de Python. Le code de dessin à introduire est disponible sur GitHub, donc si vous êtes intéressé, veuillez le télécharger.

introduction

Les articles suivants ont été trouvés liés au dessin et à la gestion des données météorologiques.

  1. Visualisez GRIB2 sur une carte avec Cartopy
  2. Visualisez grib2 sur une carte avec python (matplotlib)

C'est très facile à comprendre, mais les deux ne sont qu'une introduction au dessin d'un élément.

Je voudrais me différencier de l'article ci-dessus en introduisant le code de l'organisation des données de la figure avec un certain degré de perfection au dessin comme indiqué ci-dessous.

Cliquez ici pour que le dessin terminé soit créé cette fois. La cible capture le cas du typhon n ° 10 en 2016.

· Fig. 1 ・ Fig. 2
GPV_elem_201608210600-201608310000.gif typhoon_strength2000-2019_track_1610.png

Présentation de chaque figure

Je vais vous présenter ma propre interprétation de la figure. (... la demande peut être faible)

Nous dessinons une combinaison d'éléments météorologiques de 06h00 le 21 août 2016 à 00h00 le 31 août 2016. Le contour représente le champ d'altitude, le vecteur bleu représente le champ de vitesse du vent de 10 m / s ou plus, et la ligne diagonale + partie jaune représente la partie basse pression particulièrement forte dans le champ d'altitude spécifié.

→ Vous pouvez voir comment la puissance augmente progressivement avant d'atterrir.

Hexagrid exprime quels points les typhons des 20 dernières années (2000-2019) sont passés. De plus, le graphique montre une trajectoire de typhon spécifique, et cette fois la trajectoire du typhon n ° 10 en 2016 est tracée selon la Fig.1.

→ On peut lire à partir de l'hexagride que la plupart des typhons jusqu'à présent sont passés près du côté est de Taiwan et se sont déplacés vers le nord et ont approché le Japon, mais il a été confirmé que le typhon n ° 10 en 2016 était une route spéciale. Je peux le faire.

Données à préparer

Il peut être obtenu auprès de Institut de recherche pour une humanosphère durable (RISH). Je vais. Veuillez consulter le Centre d'assistance de la ligne météorologique pour des informations détaillées sur la grille et le délai de livraison.

Nous fournissons des données aux établissements d'enseignement et de recherche. Si vous avez fréquemment besoin de données pour les activités de l'entreprise, veuillez acheter les données directement auprès du Centre d'assistance météorologique aux entreprises et coopérer à la maintenance et au développement de l'ensemble du système de fourniture de données. (Cité sur le site)

Les données utilisées cette fois sont les données obtenues à partir du site ci-dessus, et les données de vent est-ouest, vent nord-sud, vent vertical, température et champ d'altitude sont découpées au format binaire à l'aide de wgirb2.

Le format des données découpées est décrit dans Fortran90. ... pour référence seulement.

format.f90


 integer, parameter :: nx = 720, ny = 361, np = 12
 integer, parameter :: elem = 5
 real(4)            :: gpv(elem, hgt, nx, ny)
 real(4)            :: pr(hgt)

 character(4)       :: elem_list(elem)

 data pr / 1000., 925., 850., 700., 600., 500., 400., 300., 250., 200., 150., 100. /
 data elem_list / "U   ", "V   ", "W   ", "T   ", "Z   " /

 open(*, fille=*, form='unformatted', access='direct', recl=4*nx*ny*np*elem) 
 close(*)

Il peut être obtenu par année auprès de Agence météorologique. Les informations depuis le moment du typhon jusqu'à la dépression tempérée sont décrites toutes les 6 heures.

Code d'exécution

J'espère que vous avez lu les deux codes basés sur main_driver. Si vous souhaitez l'exécuter, veuillez modifier le paramètre PATH.

Code de dessin des éléments météorologiques (Fig.1)

Tout d'abord, je présenterai le code qui dessine les éléments météorologiques.

weather_bin2plot.py


# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""

import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap

class mapping:
	
	def __init__(self):
		pass

	def gpv_data_coef(self):
		nx, ny, hgt = 720, 361, 12
		elem   = 5 
		return nx, ny, hgt, elem

	def prj_coef(self, nx, ny):
		dx, dy = 0.5, 0.5
		lon, lat = [], []
		for ix in range(nx):
			lon += [ float('{:.2f}'.format(dx*ix)) ]
		for iy in range(ny):
			lat += [ float('{:.2f}'.format(-90. + dy*iy)) ]

		X, Y = np.meshgrid(lon, lat)

		return X, Y

	def preparating_data(self, yyyy, mm, dd, hh):
		time_list = []
		num_time_list = len(time_list)
		
		month_thirtyone = [ 1, 3, 5, 7, 8, 10, 12 ]
		month_thirty    = [ 4, 6, 9, 11 ]
		month_twntynine = [ 2 ]

		while num_time_list < 40:
			time_list.append(str(yyyy) + str('%02d' % mm) + str('%02d' % dd) + str('%02d' % hh) + '00')
			hh = hh - 6
	
			if hh < 0 and dd == 1:
				mm, hh = mm - 1, 18
				if mm in month_thirty:
					dd = 30
				elif mm in month_thirtyone:
					dd = 31
				elif mm in month_twntynine:
					if yyyy % 4 == 0:
						dd = 28
					else:
						dd =29

			elif hh < 0:
				dd, hh = dd - 1, 18

			num_time_list += 1

		return time_list

	def setup_gpv_filelist(self, path, time_list):
		filelist = []
		for i_file in time_list:
			filelist.append(path + '/data')
		return filelist 

	def open_gpv_filelist(self, gpv_filelist, nx, ny, hgt, elem):
		data = [ []*i for i in range(len(gpv_filelist)) ]
		for num_gpv, gpv_file in enumerate(gpv_filelist):
			with open(gpv_file, 'rb') as ifile:
				"""
				elem: 1-u, 2-v, 3-w, 4-tmp, 5-height
				"""
				data[num_gpv] = np.fromfile(ifile, dtype='>f', sep = '').reshape(elem, hgt, ny, nx)
				print('..... Preparating data for ' + gpv_file, ', shapes :: ', data[num_gpv].shape)
		return data		

	def main_mapping_tool(self, mode, path, time_list, nx, ny, *, gpv_datalist='None', snap_step=0, level='1000'):

		fig, ax = plt.subplots()
		outpath = path + '/fig/' 
		
		lat_min, lat_max = 17, 50
		lon_min, lon_max = 120, 155
		
		mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )

		mapping.drawcoastlines(color='black', linewidth=0.5)
		mapping.drawmeridians(np.arange(0, 360, 5),  labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
		mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)

		lon_list, lat_list = self.prj_coef(nx, ny)
		x, y = mapping(lon_list, lat_list)

		if mode == 1 or mode == 2:

			if level == "1000":
				hPa, min_list = 0, [0, 200]
			elif level == "925":
				hPa, min_list = 1, [0, 600]
			elif level == "850":
				hPa, min_list = 2, [0, 750]
			elif level == "700":
				hPa, min_list = 3, [0, 1500]
			elif level == "600":
				hPa, min_list = 4, "None"
			elif level == "500":
				hPa, min_list = 5, "None"
			elif level == "400":
				hPa, min_list = 6, "None"
			elif level == "300":
				hPa, min_list = 7, "None"
			elif level == "250":
				hPa, min_list = 8, "None"
			elif level == "200":
				hPa, min_list = 9, "None"
			elif level == "150":
				hPa, min_list = 10, "None"
			elif level == "100":
				hPa, min_list = 11, "None"

			gpv_u_data = gpv_datalist[snap_step][0][hPa]
			gpv_v_data = gpv_datalist[snap_step][1][hPa]
			speed = np.sqrt(gpv_u_data*gpv_u_data + gpv_v_data*gpv_v_data)
			speed_list = list(range(0, 50, 25)) 
		
			gpv_height_data = gpv_datalist[snap_step][4][hPa]
			num_list = list(range(0, 7500, 10))

			contour = mapping.contour(x, y, gpv_height_data, linewidths=0.25, linestyles='-', levels=num_list, colors='m')
			contour.clabel(fmt='%1.1f', fontsize=6.5)

			if not min_list == "None":
				lines = mapping.contourf(x, y, gpv_height_data, min_list, alpha=0.5, hatches=['///'], lw=1., zorder=0, edgecolor='r', colors="y")

			for i_nx, i_ny in itertools.product(range(nx), range(ny)):
				if speed[i_ny][i_nx] > 10 and lon_min <= lon_list[i_ny][i_nx] <= lon_max and lat_min <= lat_list[i_ny][i_nx] <= lat_max:
					print('...... Halfway step, ', i_nx, i_ny, speed[i_ny][i_nx])
					vector = mapping.quiver(x[i_ny][i_nx], y[i_ny][i_nx], gpv_u_data[i_ny][i_nx], gpv_v_data[i_ny][i_nx], color='c', units='dots', scale=2.0, alpha=0.6)
			
			plt.title(time_list[snap_step] + ' @GSM ' + level + 'hPa' , loc='left', fontsize=10)
			plt.quiverkey(vector, 0.75, 0.9, 10, '10 m/s', labelpos='W', coordinates='figure')

			plt.savefig(outpath + 'GPV_elem_' + time_list[snap_step] + '.png')
			print('...... Saving fig :: ', outpath + 'GPV_elem_' + time_list[snap_step] + '.png')

		#plt.show()

	def main_driver(self, mode, indir, time_list, level):
		nx, ny, hgt, elem = self.gpv_data_coef()

		gpv_filelist = self.setup_gpv_filelist(indir, time_list)
		gpv_datalist = self.open_gpv_filelist(gpv_filelist, nx, ny, hgt, elem)

		if mode == 1:
			snap_step = 0
			self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)
		elif mode == 2:
			for snap_step in range(len(gpv_datalist)):
				self.main_mapping_tool(mode, indir, time_list, nx, ny, gpv_datalist=gpv_datalist, level=level, snap_step=snap_step)

if __name__ == "__main__":
	mapp = mapping()

	#input dir
	indir = '/your_directory/'

	#target time
	yyyy, mm, dd, hh = 2016, 8, 31, 0
	time_list = mapp.preparating_data(yyyy, mm, dd, hh)

	"""
	Target altitude list
	1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 150, 100 hPa
	"""
	level = '925' 

	#choose mode, if you append new func. more anynum. 
	"""
	2019.12.11
	mode 1: Normal weather element info at GPV GSM snap shot.
	mode 2: Normal weather element info at GPV GSM gif.
	"""
	
	mode = 1

	#main_driver
	mapp.main_driver(mode, indir, time_list, level)

Les points élaborés par ce code étaient les suivants.

Code de dessin de la trajectoire du typhon (Fig.2)

Ensuite, je présenterai le degré de cours de typhon au cours des 20 dernières années et le code qui a créé le cours d'un typhon spécifique.

typhoon_csv2plot.py


# -*- coding: utf-8 -*-
"""
Created on 2019.12.11
@author: Toyo_Daichi
"""

import numpy as np
import pandas as pd
import glob
import itertools
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap

class mapping:

	def __init__(self):
		pass

	def setup_csv_filelist(self, path, *, year='*'):
		filelist = []
		filelist = glob.glob(path + 'table' + str(year) + '.csv')
		return filelist

	def open_csv_filelist(self, csv_filelist, *, typhoon_number='None'):
		csv_datalist = [ []*i for i in range(len(csv_filelist)) ]
		list_num     = list(range(11)) 
		list_option  = ( 'year', 'month', 'day', 'hour(UTC)', 'typhoon_number', 'typhoon_name', 'rank','latitude', 'longitude', 'central_pressure', 'max_wind')
		for num_file, infile in enumerate(csv_filelist):
			print('..... Preparating data for ' + str(num_file) + ' ' + str(infile))
			tmp_data = pd.read_csv(infile, usecols=list_num, skiprows=1, names=list_option, sep=',')
			if typhoon_number == "None":
				tmp_lat_list  = tmp_data['latitude'].values.tolist()
				csv_datalist[num_file].append(tmp_lat_list)
				tmp_lon_list = tmp_data['longitude'].values.tolist()
				csv_datalist[num_file].append(tmp_lon_list)
				tmp_centpre_list = tmp_data['central_pressure'].values.tolist()
				csv_datalist[num_file].append(tmp_centpre_list)

			else:
				csv_datalist = []
				specific_data = tmp_data.query('typhoon_number == %s' % typhoon_number) 
				tmp_lat_list  = specific_data['latitude'].values.tolist()
				csv_datalist.append(tmp_lat_list)
				tmp_lon_list = specific_data['longitude'].values.tolist()
				csv_datalist.append(tmp_lon_list)
				tmp_centpre_list = specific_data['central_pressure'].values.tolist()
				csv_datalist.append(tmp_centpre_list)

		return csv_datalist

	def main_mapping_tool(self, path, csv_datalist, csv_specific_datalist='None', typhoon_info='None'):
		fig, ax = plt.subplots()
		outpath = path + '/fig/' 

		lat_min, lat_max = 17, 50
		lon_min, lon_max = 120, 155

		mapping = Basemap( projection='lcc', resolution="i", lat_0=35, lon_0=140, fix_aspect=(1,1), llcrnrlat=lat_min, urcrnrlat=lat_max, llcrnrlon=lon_min, urcrnrlon=lon_max )

		mapping.drawcoastlines(color='black', linewidth=0.5)
		mapping.drawmeridians(np.arange(0, 360, 5),  labels=[False, True, False, True], fontsize='small', color='gray', linewidth=0.5)
		mapping.drawparallels(np.arange(-90, 90, 5), labels=[True, False, False, True], fontsize='small', color='gray', linewidth=0.5)

		full_lat_list  = list(map(lambda x: x[0], csv_datalist))
		full_lon_list  = list(map(lambda x: x[1], csv_datalist))
	
		full_lat_list = np.sum(full_lat_list, axis=0)
		full_lon_list = np.sum(full_lon_list, axis=0)
		
		lat_list, lon_list = [], []
		for i_num in range(len(full_lat_list)):
			if(lat_min <= full_lat_list[i_num] <= lat_max and lon_min <= full_lon_list[i_num] <= lon_max):
				lat_list.append(full_lat_list[i_num])	
				lon_list.append(full_lon_list[i_num])	

		x, y = mapping(lon_list, lat_list)
		hexbin = mapping.hexbin(np.array(x), np.array(y), gridsize=[30, 30], cmap='Blues', edgecolors='gray', mincnt=8)

		if not csv_specific_datalist == "None":
			specific_lat_list, specific_lon_list = csv_specific_datalist[0], csv_specific_datalist[1]
			specific_centpre_list = csv_specific_datalist[2]

			lat_list, lon_list, centpre_list = [], [], []
			for i_num in range(len(specific_lon_list)):
				if(lat_min <= specific_lat_list[i_num] <= lat_max and lon_min <= specific_lon_list[i_num] <= lon_max):
					lat_list.append(specific_lat_list[i_num])
					lon_list.append(specific_lon_list[i_num])	
					centpre_list.append(specific_centpre_list[i_num])	

			case_x, case_y = mapping(lon_list, lat_list)

			mapping.plot(case_x, case_y, linewidth=0.5, color='red', ls='--', marker='o', ms=2)

		cbar = plt.colorbar(hexbin)
		cbar.set_label('number', fontsize=8)
		if not csv_specific_datalist == "None":
			plt.title('Course of typhoon 2000-2019' + ', ' + 'Typhoon track: ' + str(typhoon_info[1]), loc='left', fontsize=9)
		else:
			plt.title('Course of typhoon 2000-2019', loc='left', fontsize=10)
		plt.savefig(outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
		print('..... savefig ' + outpath + 'typhoon_strength2000-2019' + '_track_' + str(typhoon_info[1]) + '.png')
		#plt.show()

	def main_driver(self, indir, *, typhoon_info='None'):
		csv_filelist = self.setup_csv_filelist(indir)
		csv_datalist = self.open_csv_filelist(csv_filelist)
		"""
		csv_datalist.shape = [year][latitude][longitude][central_pressure]
		"""

		if not typhoon_info == "None": 
			print('..... Check specific case filelist')
			csv_specific_filelist = self.setup_csv_filelist(indir, year=typhoon_info[0])
			csv_specific_datalist = self.open_csv_filelist(csv_specific_filelist, typhoon_number=typhoon_info[1])
			self.main_mapping_tool(indir, csv_datalist, csv_specific_datalist=csv_specific_datalist, typhoon_info=typhoon_info)
		else:
			self.main_mapping_tool(indir, csv_datalist)

if __name__ == "__main__":
	mapp = mapping()

	#input dir
	indir = '/your_directory/'

	"""
	If you want to write a specific typhoon route, enter the typhoon information.
	typhoon_info = [year, typhoon_number]
	"""
	typhoon_info = [2016, 1610]

	#main_driver
	mapp.main_driver(indir, typhoon_info=typhoon_info)

	print('Normal END')

Les points élaborés par ce code étaient les suivants.

Environnement d'utilisation utilisé cette fois

Je suis désolé de décrire l'environnement.

requirement.txt


basemap==1.2.0 
matplotlib==3.1.1

finalement

Le code que j'ai créé est petit, mais j'ai fait preuve d'ingéniosité à ma manière, je voudrais donc le présenter à l'avenir. Si vous apprenez à mieux écrire, je vous serais reconnaissant de bien vouloir m'apprendre dans les commentaires.

Merci d'avoir lu jusqu'au bout.

Site référencé

Recommended Posts

Introduction du code de dessin pour les figures avec un certain degré de perfection des données météorologiques
[Introduction à Python] Comment obtenir l'index des données avec l'instruction for
Peut être utilisé avec AtCoder! Une collection de techniques pour dessiner du code court en Python!
Un mémorandum de méthode souvent utilisé lors de l'analyse de données avec des pandas (pour les débutants)
Recommandation de Jupyter Notebook, un environnement de codage pour les data scientists
Un diagramme de réseau a été créé avec les données du COVID-19.
Un exemple pour dessiner des points avec PIL (Python Imaging Library).
Une collection de méthodes utilisées lors de l'agrégation de données avec des pandas
Tourner un tableau de chaînes avec une instruction for (Python3)
[Python] Créer un écran pour le code d'état HTTP 403/404/500 avec Django
Gérez le chevauchement lors du dessin d'un diagramme de dispersion avec une grande quantité de données (Matplotlib, Pandas, Datashader)