Как объединить выходные данные WRF с шейпфайлом в Python?

Я хочу отобразить переменную LU_INDEX в WRF с помощью pcolormesh и добавить береговые линии из shapefile в pcolormesh. Моя координатная система WRF — LambertConformal, а shapefile — WGS84. Я думаю, что проблемы возникли из-за несоответствия двух наборов данных. Как я могу отредактировать свой код и получить четкую фигуру?

Данные проекции карты WRF

 map_proj = 'lambert',
 ref_lat   =  35.50,
 ref_lon   = 125.00,
 truelat1  =  30.0,
 truelat2  =  60.0,
 stand_lon = 125.0,

from netCDF4 import Dataset
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from wrf import (to_np, getvar, get_cartopy, latlon_coords,smooth2d)
from cartopy.feature import NaturalEarthFeature
from cartopy.feature import ShapelyFeature
from cartopy.io.shapereader import Reader
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib.colors import ListedColormap, BoundaryNorm
import matplotlib
matplotlib.use('Agg')

# Установка пути к shapefile
shp_path = "/home/chwan/data/shp/241002/gadm41_KOR_1.shp"
base_path = "/home/chwan/UCM/WPS-4.3.1/"
pdf_pages = PdfPages(f"{base_path.split('https://stackoverflow.com/')[-2]}_LU_INDEX_contour.pdf")
dates = [19]
timeidx = 0

colors = [
    "white",      # 0
    "navy",       # 1
    "grey",      # 2
    "olivedrab",     # 3
    "orange",     # 4
    "yellow",      # 5
    "green",     # 6
    "lime",  # 7
    "limegreen",        # 8
    "white",# 9
    "tan",  # 10
    "forestgreen",      # 11
    "white",    # 12
    "white",# 13
    "darkgreen",   # 14
    "olive",  # 15
    "skyblue",     # 16
    "deepskyblue",        # 17
    "seagreen",     # 18
    "white",       # 19
    "white",       # 20
    "white",    # 21
    "white",       # 22
    "chocolate",   # 23
    "white",  # 24
    "white", # 25
    "white",      # 26
    "white",  # 27
    "white",      # 28
    "white", # 29
    "white",  # 30
    "white", # 31
    "white",       # 32
    "maroon" # 33
]

colorlist = [colors[i] if i < len(colors) else "white" for i in range(34)]
non_white_colors = [color for color in colors if color != "white"]
non_white_indices = [i for i, color in enumerate(colors) if color != "white"]
cmap = ListedColormap(non_white_colors)
norm = BoundaryNorm(np.arange(0,len(colors)+1)-0.5,len(colors))
norm = BoundaryNorm(non_white_indices + [max(non_white_indices)+1],cmap.N)

globe = ccrs.Globe(ellipse="sphere", semimajor_axis=6370000,semiminor_axis=6370000)

for day in dates:
    filename = f"{base_path}geo_em.d03.nc"
    ncfile = Dataset(filename, "r")
    lu = getvar(ncfile, "LU_INDEX")
    cart_proj = get_cartopy(lu)
    lats, lons = latlon_coords(lu)
    print(cart_proj)
    print(lats)

    fig = plt.figure(figsize=(10, 10))
    ax = plt.axes(projection=ccrs.LambertConformal(globe=globe,central_longitude=125.0,central_latitude=35.5,standard_parallels=(30,60)))
    #ax.set_extent([125.5, 127.5, 34, 35.5])

    lat_mask = (to_np(lats) >= 34) & (to_np(lats) <= 35.5)
    lon_mask = (to_np(lons) >= 125.5) & (to_np(lons) <= 127.5)
    mask = lat_mask & lon_mask

    lu_np = to_np(lu)

    shpfile = Reader(shp_path)
    shapefile_feature = ShapelyFeature(shpfile.geometries(),
                                          ccrs.LambertConformal(globe=globe,central_longitude=125.0,central_latitude=35.5,standard_parallels=(30,60)),
                                          edgecolor="black", facecolor="none", linewidth=1)
    ax.add_feature(shapefile_feature)  # Добавление shapefile

    # Добавление государственных границ
    #states = NaturalEarthFeature(category="cultural", scale="10m", name="admin_0_countries",
    #                              facecolor="none", edgecolor="black")
    #ax.add_feature(states, linewidth=1, edgecolor="black")

    def get_edges(arr):
        edges = (arr[:-1] + arr[1:]) / 2.0
        return np.concatenate([[arr[0] - (edges[0] - arr[0])], edges, [arr[-1] + (arr[-1] - edges[-1])]])

    # Вычисление границ долгот и широт
    lon_edges = get_edges(to_np(lons)[0, :])
    lat_edges = get_edges(to_np(lats)[:, 0])

    mesh = ax.pcolormesh(lon_edges, lat_edges, to_np(lu), transform = ccrs.LambertConformal(globe=globe,central_longitude=125.0,central_latitude=35.5,standard_parallels=(30,60)),cmap=cmap,norm = norm, shading = 'auto', edgecolors="face")

    # Добавление цветовой шкалы
    cbar = plt.colorbar(mesh, orientation='vertical', shrink=1, pad=0.075, aspect=15)
    cbar.ax.tick_params(labelsize=11)
    cbar.set_label("ЗЕМЕЛЬНОЕ ИСПОЛЬЗОВАНИЕ", labelpad=15, rotation=270, ha="center", va="center")

    # Сохранение в PDF
    pdf_pages.savefig(fig)
    plt.close(fig)

pdf_pages.close()