Я получаю эту ошибку:
Traceback (последний вызов был последним): файл "geraEstatistica.py", строка 148, в cnplot = map.contour (x, y, resultado, clevs, cmap = plt.cm.jet) Файл "/home/caio/anaconda3/envs/showcast2/lib/python3.7/site-packages/mpl_toolkits/basemap/init.py", строка 543, в возвращении with_transform plotfun c (self, x, y, data, * args, ** kwargs) Файл "/home/caio/anaconda3/envs/showcast2/lib/python3.7/site-packages/mpl_toolkits/basemap/init .py ", строка 3569, в контурной маске = np.logical_or (ma.getmaskarray (data), xymask) ValueError: операнды не могут быть переданы вместе с фигурами (97,153) (2,2)
Когда я запускаю следующий python скрипт:
# -*- coding: utf-8 -*-
from netCDF4 import Dataset
#from mpl_toolkits.axes_grid1.inset_locator import inset_axes # Add a child inset axes to this existing axes.
from datetime import datetime, timedelta # Library to convert julian day to dd-mm-yyyy
#from cpt_convert import loadCPT # Import the CPT convert function
#from matplotlib.colors import LinearSegmentedColormap # Linear interpolation for color maps
import matplotlib.pyplot as plt # Plotting library
#import matplotlib.colors # Matplotlib colors
import numpy as np # Scientific computing with Python
#import cartopy, cartopy.crs as ccrs # Plot maps
#import cartopy.io.shapereader as shpreader # Import shapefiles
#import time as t # Time access and conversion
import math # Import math
#from matplotlib.image import imread # Read an image from a file into an array
import os # Miscellaneous operating system interfaces
#import sys # Import the "system specific parameters and functions" module
#from cartopy.feature.nightshade import Nightshade # Draws a polygon where there is no sunlight for the given datetime.
#from html_update import update # Update the HTML animation
from remap import remap # Import the Remap function
from mpl_toolkits.basemap import Basemap
#from netCDF4 import Dataset as open_ncfile
import re
global extent_
global lat
global lon
dpi = 100
def interpreta(path): #interpreta os dados da matriz ACMF Clear Sky Mask
#start = t.time()
#nc = open_ncfile('netsacmf_jan_2020/'+str(path))
#print(nc.variables)
# Image path
path = "netsacmf_jan_2020/"+str(path)+"_SEC"
# Remove the composite id
path_proc = path[:-4]
# Read the image
file = Dataset(path_proc)
# Read the satellite
satellite = getattr(file, 'platform_ID')
# Read the main variable
variable = list(file.variables.keys())[0]
# Read the product name
prod_name = getattr(file, 'title')
# Desired resolution
resolution = 8
# Read the resolution
band_resolution_km = getattr(file, 'spatial_resolution')
band_resolution_km = float(band_resolution_km[:band_resolution_km.find("km")])
# Division factor to reduce image size
f = math.ceil(float(resolution / band_resolution_km))
#print(resolution)
#print(band_resolution_km)
#print(f)
# Read the central longitude
longitude = file.variables['goes_imager_projection'].longitude_of_projection_origin
# Read the semi major axis
a = file.variables['goes_imager_projection'].semi_major_axis
# Read the semi minor axis
b = file.variables['goes_imager_projection'].semi_minor_axis
# Calculate the image extent
h = file.variables['goes_imager_projection'].perspective_point_height
x1 = file.variables['x_image_bounds'][0] * h
x2 = file.variables['x_image_bounds'][1] * h
y1 = file.variables['y_image_bounds'][1] * h
y2 = file.variables['y_image_bounds'][0] * h
# Reading the file time and date
add_seconds = int(file.variables['time_bounds'][0])
date = datetime(2000,1,1,12) + timedelta(seconds=add_seconds)
date_formated = date.strftime('%Y-%m-%d %H:%M UTC')
date_file = date.strftime('%Y%m%d%H%M')
year = date.strftime('%Y')
month = date.strftime('%m')
day = date.strftime('%d')
hour = date.strftime('%H')
minutes = date.strftime('%M')
# Choose the visualization extent (min lon, min lat, max lon, max lat)
extent_ = [float("-54.0"), float("-25.0"), float("-43.0"), float("-18.0")]
# Call the reprojection funcion
grid = remap(path_proc, variable, extent_, resolution, h, a, b, longitude, x1, y1, x2, y2)
# Read the data returned by the function
data = grid.ReadAsArray()
# Convert from int16 to uint16
data = data.astype(np.float64)
return data
print('- Mapa de janeiro 2020 - cobertura de nuvem - \n Script iniciado! \n gerando matrizes binárias... \n calculando soma de matrizes... ')
global arquivos;
global soma;
soma = np.array([], dtype = float)
for arquivo in os.walk('netsacmf_jan_2020/'): #lê o diretório que contém os nets e transforma em listas com os nomes dos arquivos
arquivos = list(arquivo[2])
for arq in arquivos:
if not re.search('aux.xml', arq):
net_receptus = interpreta(arq) #manda para a função interpreta para realizar a leitura do nets e transformar na matriz binária net_receptus
if (len(soma) == 0):
soma = net_receptus
else:
soma = np.add(soma,net_receptus)
print("calculando estatística final de janeiro...")
qttd = len(arquivo[2])
resultado = (soma / qttd)
lat = np.array([90, -90], dtype = float)
lon = np.array([0, 357.5], dtype = float)
fig = plt.figure(figsize=(1100/dpi, 1100/dpi), dpi=dpi)
ax = fig.add_axes([0.1,0.1,0.8,0.9])
map = Basemap(projection='cyl',llcrnrlat= -90.,urcrnrlat= 90.,
resolution='c', llcrnrlon=-180.,urcrnrlon=180.)
map.drawcoastlines()
map.drawstates()
map.drawcountries()
map.drawparallels(np.arange( -90., 90.,30.),labels=[1,0,0,0],fontsize=10)
map.drawmeridians(np.arange(-180.,180.,30.),labels=[0,0,0,1],fontsize=10)
x, y = map(*np.meshgrid(lon,lat))
clevs = np.arange(0,100,1)
cnplot = map.contour(x,y,resultado,clevs,cmap=plt.cm.jet)
plt.title('Mapa de indice de combertura de nuvens - Estado de São Paulo - Janeiro 2020')
plt.show()
plt.savefig('mapa_de_indice_de_cobertura_de_nuvens.png', bbox_inches='tight', dpi=dpi)
Цель этого скрипта - построить матрицу результатов, которая представляет собой матрицу сумм, деленную на количество файлов netCDF с прочитанными данными покрытия облаков. Матричная сумма - это сумма всех двоичных матриц, которые генерируются файлами netCDF.
[[0,1,0,1,1,0]] + [[0,1,0,0,1,1]] + ....... = [[20,40,0,20,60,10]] / amount of binary arrays or amount of netCDF's files = matrix result with percentages.
Я пытаюсь отобразить этот результат матрицы с процентами на карте, которая выглядит как на рисунке ниже (только то, что я делаю для карты состояний SP):
image
для этого я использую функции countr() и countrf() в matplotlib. но я все еще получаю эту трассировку.
Я использую файлы netCDF с ACMF с маской облачного покрова, доступной через amazon AWS.
Существует несколько таких файлов:
OR_ABI-L2-ACMF-M6_G16_s202000100216_e202000109524_c202000110376.nc каждый генерирует двоичный массив. 0 = нет облаков и 1 = с облаками.
вот папка с файлами netCDFs, поэтому вы можете запустить скрипт:
ссылка google drive