Fonctions utiles

mail@pastecode.io avatar
unknown
python
2 years ago
24 kB
6
Indexable
Never
import numpy as np
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset as NetCDFFile
import matplotlib.pyplot as plt
from matplotlib import ticker
import skfuzzy as fuzz
from skfuzzy import control as ctrl
import imageio
import os
from geopy import distance
from scipy.optimize import curve_fit
import scipy
import similaritymeasures
from sewar.full_ref import mse, rmse, psnr, uqi, ssim, ergas, scc, rase, sam, msssim, vifp
import math


# Pour pouvoir extraire les données:
# longitude / latitude / time / d = divergence / r = humidity / t = temperature / u=uwind / v=vwind / vo=vortcity
# sp = surface_pressure

#Pour déssiner une map pour une variable:
def plotNC(path, var, unite, titre,jour,heure):
    nc = NetCDFFile(path) # Ouverture de la BDD

    #Extraction des tableaux de données:
    lat = nc.variables['latitude'][:]
    lon = nc.variables['longitude'][:]
    time = nc.variables['time'][:]
    vorticity = nc.variables['vo'][:]
    var = nc.variables[var][:] #Ou une autre variable au choix dans la définition de la fonction

    nc.close()

    #Création de la map
    map = Basemap(width=5000000, height=3500000,
                  resolution='l', projection='cyl',
                  llcrnrlon=lon.min(), llcrnrlat=lat.min(), urcrnrlon=lon.max(), urcrnrlat=lat.max(), lat_0=lat.mean(),
                  lon_0=lon.mean())
    lons, lats = np.meshgrid(lon, lat)
    xi, yi = map(lons, lats)

    map.drawmapboundary(fill_color='aqua')
    map.fillcontinents(color='coral', lake_color='aqua')
    map.drawcoastlines()

    parallels = np.arange(lat.min(), lat.max(), 5.)  # make latitude lines every 5 degrees from xx to xx
    meridians = np.arange(lon.min(), lat.min(), 5.)  # make longitude lines every 5 degrees from xx to xx
    map.drawparallels(parallels, labels=[1, 0, 0, 0], fontsize=5)
    map.drawmeridians(meridians, labels=[0, 0, 0, 1], fontsize=5)


    date_exact = (jour-1) * 24 + heure #Date exacte, décalage de -1 jour car la basse de données commence au jour 0

    #Création des données en couleur
    cs = map.pcolor(xi, yi, np.squeeze(var[date_exact, :, :]),vmin = var.min(), vmax = var.max())
    cbar = map.colorbar(cs, location='bottom', pad=0.2,)
    if titre == 'Divergence':
        tick_locator = ticker.MaxNLocator(nbins=5)
        cbar.locator = tick_locator
        cbar.update_ticks()
    cbar.set_label(unite,size=7)
    plt.title(titre,fontsize = 10) #Pour déssiner les map

#Pour déssiner les map pour toutes les variables, avec export en photos:
def plotTC(path):

    nc = NetCDFFile(path)
    j = [23,24,25,26,27,28,29,30,31]
    h = [10,23]

    vorticity = np.array(nc.variables['vo'][:])
    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])

    # np.shape(vorticity) = (744,93,117)
    # np.shape(latitude) = (93,)
    # np.shape(longitude) = (117,)
    # L'indice sera la latitude et la longitude. Cependant,ce ne sont pas les vraies lat et lon.
    # Il faut prendre la matrice nc.variable de la latitude puis l'évaluer en l'index trouvé pour avoir la vraie latitude.
    # Et pareil pour la longitude.

    time = np.shape(vorticity)[0]

    x = []
    y = []

    for i in range(time):
        a = np.unravel_index(np.argmax(vorticity[i,:,:]), vorticity[i,:,:].shape) #prend le max de la vorticité qui correspond à la position du cyclone
        latitude_coord, longitude_coord = a[0],a[1]
        x.append(latitude[latitude_coord])
        y.append(longitude[longitude_coord])

    #Et pour chaque jour on crée les images:

    i = 0

    for jour in j:
        for heure in h:
            plt.rcParams.update({'font.size': 5})
            plt.subplot(2, 3, 1)
            plotNC(path, 'r', '%', 'Humidité relative',jour,heure)
            plt.subplot(2, 3, 2)
            plotNC(path, 'd', 'jsp', 'Divergence',jour,heure)
            plt.subplot(2, 3, 3)
            plotNC(path, 'vo', 'jsp', 'Vorticité',jour,heure)
            plt.subplot(2, 3, 4)
            plotNC(path, 't', 'Kelvin', 'Température',jour,heure)
            plt.subplot(2, 3, 5)
            plotNC(path, 'u', 'm/s', 'Uwind',jour,heure)
            plt.subplot(2, 3, 6)
            plotNC(path, 'v', 'm/s', 'Vwind',jour,heure)
            plt.savefig(f"{i}.png", dpi=300)
            plt.clf()
            i = i+1

    #Pour faire une belle séquence des images:
    graphs = []
    for i in range(0,16):
        graphs.append(imageio.imread(f"{i}.png"))

    imageio.mimsave(f"{os. getcwd()}.gif",graphs)

    plt.show()

#Nous donne en sortie la position du TC à un moment donné (Latitude, Longitude)

def position_TC(jour_TC, heure_TC, path):
    nc_TC = NetCDFFile(path)
    date = (jour_TC-1) * 24 + heure_TC
    v = np.array(nc_TC.variables['v'][:])
    u = np.array(nc_TC.variables['u'][:])
    vor = np.array(nc_TC.variables['vo'][:])
    lat = np.array(nc_TC.variables['latitude'][:])
    lon = np.array(nc_TC.variables['longitude'][:])

    u = u[date, :, :]
    v = v[date, :, :]
    vor = vor[date, :, :]

    x = scipy.integrate.cumtrapz(v, lon, axis=1, initial=0)
    y = scipy.integrate.cumtrapz(u, lat, axis=0, initial=0)
    psi_vor = (y - x)
    psi_vor = psi_vor / np.min(psi_vor) + vor / np.min(vor)

    latitude = np.array(nc_TC.variables['latitude'][:])
    longitude = np.array(nc_TC.variables['longitude'][:])
    # On trouve la vorticité max à cette date, dans une certaine zone
    a = np.unravel_index(np.nanargmin(psi_vor), psi_vor.shape)

    # Puis on envoie les coordonnées
    # print(f"(Lon,Lat) = {latitude[a[0]], longitude[a[1]]}")

    nc_TC.close()

    return latitude[a[0]], longitude[a[1]]

def ffposition_TC(jour,heure,path):
    nc = NetCDFFile(path)
    date_exact = (jour-1) * 24 + heure
    vwind = np.array(nc.variables['v'][:])
    uwind = np.array(nc.variables['u'][:])
    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])
    vo = np.array(nc.variables['vo'][:])
    x = scipy.integrate.cumtrapz(vwind[date_exact, :, :], longitude, axis=1, initial=0)
    y = scipy.integrate.cumtrapz(uwind[date_exact, :, :], latitude, axis=0, initial=0)
    psi_vor = (y - x)
    psi_vor = psi_vor / np.min(psi_vor) + vo / np.min(vo)

    a = np.unravel_index(np.nanargmin(psi_vor), psi_vor.shape)

    # Puis on envoie les coordonnées
    return latitude[a[0]], longitude[a[1]]

def fposition_TC(jour, heure, path):
    nc = NetCDFFile(path)
    date_exact = (jour - 1) * 24 + heure

    vorticity = np.array(nc.variables['vo'][:])
    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])

    # On trouve la vorticité max à cette date
    a = np.unravel_index(np.argmax(vorticity[date_exact, :, :]), vorticity[date_exact, :, :].shape)

    # Puis on envoie les coordonnées
    return (latitude[a[0]], longitude[a[1]])  # Return Lat,Lon

#Nous donne en sortie le paramètre de coriolis
def parametre_coriolis(latitude):

    latitude = (latitude*np.pi)/180

    return 2*0.72921*(10**(-4))*np.sin(latitude)

#Nous donne en sortie la taille du Cyclone (diamètre le plus grand de l'ellipse)
def taille_TC(jour, heure, seuil, path):
    nc = NetCDFFile(path)
    date_exact = (jour - 1) * 24 + heure

    # np.shape(vorticity) = (744,93,117)
    # np.shape(latitude) = (93,)
    # np.shape(longitude) = (117,)

    vorticity = np.array(nc.variables['vo'][:])
    latitude = np.array(nc.variables['latitude'][:])
    longitude = np.array(nc.variables['longitude'][:])

    # print(longitude[80]) gives -74.0

    # Limit longitude study up to -74.0 to the right to ignore high vorticity where there is no cyclone; can be observed on the .gif

    vo_max = np.unravel_index(np.argmax(vorticity[date_exact, :, 0:74]), vorticity[date_exact, :, 0:74].shape)
    lon_TC = vo_max[1]
    lat_TC = vo_max[0]

    s = seuil
    vorticity[vorticity < s] = 0

    # Calculating min and max lon and lat

    # Hovering the latitude, fixing longitude of current position

    lat_vo_min = 0  # Initialising values
    lat_vo_max = 0  # Initialising values
    for i in range(93):
        if vorticity[
            date_exact, i, lon_TC] > 0:  # Taking the first non zero value index as the start point of the TC latitude
            lat_vo_min = i
            break

    if lat_vo_min == 0:  # If lat_vo_min==0 it means that all values were below s
        lat_vo_max = 0
    else:
        for i in range(lat_vo_min, 93):
            if vorticity[
                date_exact, i, lon_TC] == 0:  # Taking the last non zero value index as the end point of the TC latitude
                lat_vo_max = i - 1
                break

    # Hovering the longitude, fixing latitude of current position

    lon_vo_min = 0  # Initialising values
    lon_vo_max = 0  # Initialising values
    for i in range(117):
        if vorticity[
            date_exact, lat_TC, i] > 0:  # Taking the first non zero valueindex as the start point of the TC longitude
            lon_vo_min = i
            break

    if lon_vo_min == 0:
        lon_vo_max = 0
    else:
        for i in range(lon_vo_min, 93):
            if vorticity[
                date_exact, lat_TC, i] == 0:  # Taking the last non zero value index as the end point of the TC longitude
                lon_vo_max = i - 1
                break

    # print(f"LAT : Min:{lat_vo_min},Max:{lat_vo_max}")
    # print(f"LON : Min{lon_vo_min},Max:{lon_vo_max}")

    # ---------------------------------------------------------------------------------

    lat1, lon1, lat2, lon2, R = lat_vo_min, lon_vo_min, lat_vo_max, lon_vo_max, 6373.0

    # Distance between lon_vo_min and lon_vo_max, lat_TC constant
    coordinates_from = [lat_TC, lon1]
    coordinates_to = [lat_TC, lon2]

    distance_geopy = distance.distance(coordinates_from, coordinates_to).km
    distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km

    longueur_longitudinale = np.rint((distance_geopy + distance_geopy_great_circle) / 2)
    # print('Longueur longitudinale', longueur_longitudinale)

    # Distance between lat_vo_min and lat_vo_max, lon_TC constant

    coordinates_from = [lat1, lon_TC]
    coordinates_to = [lat2, lon_TC]

    distance_geopy = distance.distance(coordinates_from, coordinates_to).km
    distance_geopy_great_circle = distance.great_circle(coordinates_from, coordinates_to).km

    longueur_latitudinale = np.rint((distance_geopy + distance_geopy_great_circle) / 2)

    # print('Longueur latitudinale', longueur_latitudinale)

    # print('AIRE', longueur_latitudinale*longueur_longitudinale,'M2')

    res = max(longueur_latitudinale, longueur_longitudinale)
    return res

#Calcule le module de la vitesse en tout point à tout moment
def calcul_pth(path):
    nc = NetCDFFile(path)
    uwind = np.array(nc.variables['u'][:])
    vwind = np.array(nc.variables['v'][:])
    pth = np.sqrt(np.power(uwind, 2) + np.power(vwind, 2))

    return(pth)

#Donne en sortie la vitesse du TC (utilise la matrice créée juste avant avec calcul_pth)
def vitesse_TC(jour,heure,pth):

    res = np.max(pth[(jour-1)*24+heure])

    return(res)  # Pythagore

#Donne la sortie la variation de vitesse du TC, à distance d'une heure)
def variation_vitesse_TC(jour,heure,pth):

    rep = vitesse_TC(jour,heure,pth)-vitesse_TC(jour,heure-1,pth)

    return(rep)

#Donne en sortie la pression à la surface où se trouve le TC, à un moment donné
def pression_TC(jour,heure,path,pressure_path):

    nc = NetCDFFile(path)
    pressure = np.array(NetCDFFile(pressure_path).variables['sp'][:])
    date_exact = (jour-1) * 24 + heure
    vorticity = np.array(nc.variables['vo'][:])


    # Puis on envoie les coordonnées
    lat,lon = position_TC(jour,heure,path)

    return(np.rint(pressure[date_exact,lat,lon]))


def variation_pression_TC(jour,heure,path,pressure_path):
    return pression_TC(jour,heure,path,pressure_path)-pression_TC(jour,heure-1,path,pressure_path)


def relation_pression_taille_TC(path,pression):

    t = []
    taille = []
    pressure = []

    for i in range(23, 32):
        for h in range(24):
            t.append((i - 1) * 24 + h)
            pressure.append(pression_TC(i, h, path, pression))
            taille.append(taille_TC(i, h, 0.0006, path))

    fig, ax1 = plt.subplots()
    color = 'tab:red'
    ax1.set_xlabel('Temps')
    ax1.set_ylabel('Pression (Pa)', color=color)
    ax1.plot(t, pressure, color=color)
    ax1.tick_params(axis='y', labelcolor=color)

    ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

    color = 'tab:blue'
    ax2.set_ylabel('Taille (Km)', color=color)  # we already handled the x-label with ax1
    ax2.plot(t, taille, color=color)
    ax2.tick_params(axis='y', labelcolor=color)

    fig.tight_layout()  # otherwise the right y-label is slightly clipped
    plt.show()


def relation_pression_vitesse_TC(path,pression):
    pth = calcul_pth(path)
    t = []
    taille = []
    pressure = []

    for i in range(23, 32):
        for h in range(24):
            t.append((i - 1) * 24 + h)
            pressure.append(pression_TC(i, h, path, pression))
            taille.append(vitesse_TC(i, h, pth))

    fig, ax1 = plt.subplots()
    color = 'tab:red'
    ax1.set_xlabel('Temps')
    ax1.set_ylabel('Pression (Pa)', color=color)
    ax1.plot(t, pressure, color=color)
    ax1.tick_params(axis='y', labelcolor=color)

    ax2 = ax1.twinx()  # instantiate a second axes that shares the same x-axis

    color = 'tab:blue'
    ax2.set_ylabel('Vitesse (Km/h)', color=color)  # we already handled the x-label with ax1
    ax2.plot(t, taille, color=color)
    ax2.tick_params(axis='y', labelcolor=color)

    fig.tight_layout()  # otherwise the right y-label is slightly clipped
    plt.show()


def regression_trajectoire_TC(debut, fin, path):
    x_full = []
    y_full = []

    # Données full
    for i in range(23, 32):
        for h in range(24):
            pos1 = position_TC(i, h, path)[0]
            pos2 = position_TC(i, h, path)[1]
            x_full.append(pos1)
            y_full.append(pos2)

    # Données entrainement
    x = []
    y = []

    for i in range(debut, fin+1):
        for h in range(24):
            pos1 = position_TC(i, h, path)[0]
            pos2 = position_TC(i, h, path)[1]
            x.append(pos1)
            y.append(pos2)

    plt.scatter(y, x)
    plt.show()

    def reg_poly(x, y, x_full, y_full):

        def objective(x, a, b, c):
            return a * (x ** 2) + b * x + c

        params = np.polyfit(x,y,2)
        a = params[0]
        b = params[1]
        c = params[2]

        y_test = []
        for pt in x_full:
            y_test.append(objective(pt, a, b, c))

        y_test = np.array(y_test)

        plt.scatter(y_test, x_full)
        plt.show()

        erreur = []
        for i in range(len(y_full)):
            e = (y_full[i] - y_test[i]) ** 2
            erreur.append(np.sqrt(e))

        print(f"Erreur poly en moyenne:{np.sum(erreur)}")
    reg_poly(x, y, x_full, y_full)

    def reg_lin(x, y, x_full):
        def objective(x, a, b):
            return a * x + b

        params = curve_fit(objective, x, y)
        a, b = params[0]

        y_test = []
        for pt in x_full:
            y_test.append(objective(pt, a, b))

        plt.scatter(y_test, x_full)
        plt.show()

        erreur = []
        for i in range(len(y_full)):
            e = (y_full[i] - y_test[i]) ** 2
            erreur.append(np.sqrt(e))

        print(f"Erreur lin en moyenne:{np.sum(erreur)}")
    reg_lin(x, y, x_full)


def distance_minimale(jour,path):

    x_full = []
    y_full = []

    # Données full
    for i in range(23, 32):
        for h in range(24):
            pos1 = position_TC(i, h, path)[0]
            pos2 = position_TC(i, h, path)[1]
            x_full.append(pos1)
            y_full.append(pos2)

    # Données entrainement
    x = []
    y = []

    for i in range(23, jour):
        for h in range(24):
            pos1 = position_TC(i, h, path)[0]
            pos2 = position_TC(i, h, path)[1]
            x.append(pos1)
            y.append(pos2)

    plt.scatter(y,x,label='Entrainement')
    plt.legend()
    plt.show()

    def reg_poly(x, y, x_full, y_full):

        def objective(x, a, b, c):
            return a * (x ** 2) + b * x + c

        params = np.polyfit(x,y,2)
        a = params[0]
        b = params[1]
        c = params[2]

        y_test = []
        for pt in x_full:
            y_test.append(objective(pt-5, a, b, c)-6)

        y_test = np.array(y_test)

        plt.scatter(y_test, x_full, label = "Estimation poly2")
        plt.scatter(y_full, x_full, label="Real")
        plt.legend()
        plt.show()

        erreur = []
        for i in range(len(y_full)):
            e = (y_full[i] - y_test[i]) ** 2
            erreur.append(np.sqrt(e))

        print(f"Erreur poly totale:{np.sum(erreur)}")

        xy_entr_futur=[]

        for i in range(len(y_full)):
            if x_full[i]>position_TC(jour,0,path)[0]:
                pos1 = x_full[i]
                pos2 = objective(pos1-5,a,b,c)-6
                xy_entr_futur.append((pos1,pos2))

        return xy_entr_futur

    def reg_lin(x, y, x_full,y_full):
        def objective(x, a, b):
            return a * x + b

        params = curve_fit(objective, x, y)
        a, b = params[0]

        y_test = []
        for pt in x_full:
            y_test.append(objective(pt, -1*a, b)+36)

        y_test = np.array(y_test)
        plt.scatter(y_test, x_full, label = "Estimation lin")
        plt.scatter(y_full, x_full, label="Real")
        plt.legend()
        plt.show()

        erreur = []
        for i in range(len(y_full)):
            e = (y_full[i] - y_test[i]) ** 2
            erreur.append(np.sqrt(e))

        xy_entr_futur = []

        for i in range(len(y_full)):
            if x_full[i]>position_TC(jour,0,path)[0]:
                pos1 = x_full[i]
                pos2 = objective(pos1,-1*a,b)+36
                xy_entr_futur.append((pos1,pos2))

        print(f"Erreur lin totale:{np.sum(erreur)}")

        return xy_entr_futur

    xy_entr_poly = reg_poly(x, y, x_full, y_full)
    xy_entr_lin = reg_lin(x,y,x_full,y_full)

    #Soit pos_current la position à l'instant où on arrête les données d'entrainement
    #Et à partir duquel on souhaite faire les prédictions
    pos_current = position_TC(jour,0,path) #(Lat,Lon)

    # ENTRAINEMENT POLY2

    temp = []
    for i in range(len(xy_entr_poly)):
        temp.append(pos_current)

    distance_tous_poly2=[]
    for i in range(len(xy_entr_poly)):
        distance_tous_poly2.append(math.dist(xy_entr_poly[i],temp[i]))

    distance_min_coord_poly2 = distance_tous_poly2.index(min(distance_tous_poly2))

    print('Poly2:')
    print(f"Pt de distance min coord: {distance_min_coord_poly2}")
    distance_min = xy_entr_poly[distance_min_coord_poly2]
    print(f"Pt de distance min: {distance_min}")

    # ENTRAINEMENT LIN

    temp = []
    for i in range(len(xy_entr_lin)):
        temp.append(pos_current)

    distance_tous_lin = []
    for i in range(len(xy_entr_lin)):
        distance_tous_lin.append(math.dist(xy_entr_lin[i], temp[i]))

    distance_min_coord_lin = distance_tous_lin.index(min(distance_tous_lin))

    print('Lin:')
    print(f"Pt de distance min coord: {distance_min_coord_lin}")
    distance_min = xy_entr_lin[distance_min_coord_lin]
    print(f"Pt de distance min: {distance_min}")


def plot_array(array, lon_plot, lat_plot):


    map_plot = Basemap(width=5000000, height=3500000,
                       resolution='l', projection='cyl',
                       llcrnrlon=lon_plot.min(), llcrnrlat=lat_plot.min(), urcrnrlon=lon_plot.max(),
                       urcrnrlat=lat_plot.max(), lat_0=lat_plot.mean(),
                       lon_0=lon_plot.mean())
    lons, lats = np.meshgrid(lon_plot, lat_plot)
    xi, yi = map_plot(lons, lats)

    map_plot.drawmapboundary(fill_color='aqua')
    map_plot.fillcontinents(color='coral', lake_color='aqua')
    map_plot.drawcoastlines()
    cs = map_plot.pcolor(xi, yi, array)
    cbar = map_plot.colorbar(cs, location='bottom', pad="10%")


def met_to_deg(x):
    return x / 111000


def deg_to_met(x):
    return x * 111000


def cart_to_polar(x0, y0, x1, y1):
    # Convertit des coordonnées cartésiennes (x,y) en polaires (r,theta), dans le repère (o,x0,y0)
    # Cette fonction fonctionne aussi avec des matrices
    x0 = met_to_deg(x0)
    y0 = met_to_deg(y0)
    x1 = met_to_deg(x1)
    y1 = met_to_deg(y1)
    x_ref = x1 - x0
    y_ref = y1 - y0
    x_ref, y_ref = np.meshgrid(x_ref, y_ref)
    r = np.sqrt(deg_to_met(x_ref) ** 2 + deg_to_met(y_ref) ** 2)
    angle = np.arctan2(x_ref, y_ref)

    return r, angle


def polar_to_cart(r,angle):
    angle = np.radians(angle)
    x = r*np.cos(angle)
    y = r*np.sin(angle)
    return x,y


def angle(jour,heure,data_path):
    # print(position_TC(jour, heure-12, data_path))
    # print(position_TC(jour, heure, data_path))
    nc = NetCDFFile(data_path)
    date_exact = (jour - 1) * 24 + heure
    u = np.array(nc.variables['u'][:])
    lon = np.array(nc.variables['longitude'][:])
    lat = np.array(nc.variables['latitude'][:])

    y0, x0 = position_TC(jour, heure-6, data_path)
    y1, x1 = position_TC(jour, heure, data_path)
    y0 = deg_to_met(y0)
    x0 = deg_to_met(x0)
    y1 = deg_to_met(y1)
    x1 = deg_to_met(x1)
    r, angle = cart_to_polar(x0, y0, x1, y1)
    angle = np.degrees(angle[0][0])

    # plot_array(u[date_exact,:,:],lon,lat)
    # plt.plot(met_to_deg(x0), met_to_deg(y0), ms=5, marker='o', markeredgecolor='red')
    # plt.plot(met_to_deg(x1), met_to_deg(y1), ms=8, marker='x', markeredgecolor='green')
    # plt.show()

    return np.round(angle, decimals = 1)


def entree_SF4(jour,heure,data_path):

    return angle(jour,heure,data_path)-angle(jour,heure-6,data_path)


def xy_to_indx_lonlat(x, y, lon_map, lat_map):
    # On prend l'index ou (x,y)=(lon,lat), avec une tolérance de 25 km,
    # sachant que l'écart entre 2 valeurs de (lon,lat) est d'env 25 km aussi
    idx_lon = np.where(np.isclose(lon_map, x, atol=0.5))[-1][-1]  # sachant que la résolution est de env 25000 km
    idx_lat = np.where(np.isclose(lat_map, y, atol=0.5))[-1][-1]
    return idx_lon, idx_lat


def output_fuzz(csq,output_var):
    output = np.zeros(1)
    first_rule = True
    for membi in output_var.terms:
        csq_i = '['+str(output_var.terms[membi])+']'
        if csq[csq_i] != 0:
            x_csq = fuzz.fuzzymath.interp_universe(output_var.universe,output_var.terms[membi].mf , 1)
            if first_rule:
                output = np.array([[np.max(x_csq),csq[csq_i]]])

                first_rule = False
            else:
                output = np.append(output,[[(np.min(x_csq)),csq[csq_i]]],axis=0)
    return output