Monday, 13 June 2022

Advance a Basemap Plot

Hi I would like to interpolate a map and do this with the package Pykrige. For this I create the Meshgrid:

def get_data(df):
    return {
        "lons": df['Longitude'],
        "lats": df['Latitude'],
        "values": df['O18'],
        "alts": df['Altitude'],
    }

def extend_data(data):
    return {
        "lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
        "lats":  np.concatenate([data["lats"], data["lats"], data["lats"]]),
        "values":  np.concatenate([data["values"], data["values"], data["values"]]),
        "alts": np.concatenate([data["alts"], data["alts"], data["alts"]]),
    }

def generate_grid(data, basemap, delta=1):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta) # dont extrapolate towards the poles
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

The interpolation then looks like this:

def interpolate(data, grid):
    UK = UniversalKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        specified_drift = data["alts"], 
        # nlags=100, 
    )
    return UK.execute("grid", grid["lon"], grid["lat"])

With the command I plot the interpolated data in a map and use contourf.

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(grid["x"], grid["y"],  interpolation,100, cmap='jet', ) #plot the data on the map. plt.cm.RdYlBu_r
    color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%") 

However, the result for Europe is very disappointing. (White Dots: Messuring Stations) enter image description here How can I improve the code so that the representation looks even more detailed? Just as a note, the colour range is already very fine-tuned and unfortunately cannot be used for this.

As you can see, there are many stations in Europe, all with different values. Therefore, the Contourf plot should be much more accurate. For example, I can give this example: https://images.app.goo.gl/sPTMM6xe5Wj5QtGN8

my code is:

from traceback import print_tb
import numpy as np
from pykrige.ok import OrdinaryKriging
from pykrige.kriging_tools import write_asc_grid
import pykrige.kriging_tools as kt
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.patches import Path, PathPatch
import pandas as pd

def load_data():
    df = pd.read_csv(r"C:/Users/Oliver Weisser/Desktop/Bachelor/Programm/Daten/Classified/Mean_O18/euras_2000.csv")
    return(df)

def get_data(df):
    return {
        "lons": df['Longitude'],
        "lats": df['Latitude'],
        "values": df['O18a'],
    }

def extend_data(data):
    # Copy data to the "left" and "right" to allow for interpolation to the edges of the map
    return {
        "lons": np.concatenate([np.array([lon-360 for lon in data["lons"]]), data["lons"], np.array([lon+360 for lon in data["lons"]])]),
        "lats":  np.concatenate([data["lats"], data["lats"], data["lats"]]),
        "values":  np.concatenate([data["values"], data["values"], data["values"]]),
    }

def generate_grid(data, basemap, delta=1):
    grid = {
        'lon': np.arange(-180, 180, delta),
        'lat': np.arange(np.amin(data["lats"]), np.amax(data["lats"]), delta) # dont extrapolate towards the poles
    }
    grid["x"], grid["y"] = np.meshgrid(grid["lon"], grid["lat"])
    grid["x"], grid["y"] = basemap(grid["x"], grid["y"])
    return grid

def interpolate(data, grid):
    OK =OrdinaryKriging(
        data["lons"],
        data["lats"],
        data["values"],
        variogram_model='exponential',
        nlags=40, 
        pseudo_inv=True,
        weight=True,
        verbose=True,
    )
    return OK.execute("grid", grid["lon"], grid["lat"])

def prepare_map_plot():
    figure, axes = plt.subplots(figsize=(10,10))
    basemap = Basemap(projection='robin', lon_0=0, lat_0=0, resolution='h',area_thresh=1000,ax=axes) 
    basemap.drawcoastlines() 
    basemap.drawparallels(np.arange(-90.,120.,30.))
    basemap.drawmeridians(np.arange(0.,420.,60.))
    return figure, axes, basemap

def plot_mesh_data(interpolation, grid, basemap):
    colormesh = basemap.contourf(grid["x"], grid["y"],  interpolation,32, cmap='RdBu_r', ) #plot the data on the map. plt.cm.RdYlBu_r
    color_bar = basemap.colorbar(colormesh,location='bottom',pad="10%") 

df = load_data()
base_data = get_data(df)
figure, axes, basemap = prepare_map_plot()
grid = generate_grid(base_data, basemap, 1)
extended_data = extend_data(base_data)
interpolation, interpolation_error = interpolate(extended_data, grid)
plot_mesh_data(interpolation, grid,basemap)
plt.show()



from Advance a Basemap Plot

No comments:

Post a Comment