import cartopy.crs as ccrs
import pyproj
from xgrads import CtlDescriptor, open_CtlDataset
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from shapely.geometry import Point  # 经纬度转换为点
from xgrads import CtlDescriptor, open_CtlDataset, get_data_projection

# load data
ctl_path = 'grid.d1.ctl'
ctl = CtlDescriptor(file=ctl_path)
da = open_CtlDataset(ctl_path)
da = da.isel(time=0)


# central_latitude
central_latitude = 30.1

globe = ccrs.Globe(
    ellipse='sphere', semimajor_axis=6370000, semiminor_axis=6370000)
pdef = ctl.pdef

wgs_proj = pyproj.Proj(proj='latlong', datum='WGS84')

temp = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat)
)

e, n = pyproj.transform(wgs_proj, temp, pdef.lonref, pdef.latref)

crs = get_data_projection(ctl, globe=globe)

crs_new = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat),
    false_easting=(pdef.isize - 1) / 2. * pdef.dx - e,
    false_northing=(pdef.jsize - 1) / 2. * pdef.dy - n
)

crs_new2 = ccrs.LambertConformal(
    globe=globe,
    central_latitude=central_latitude,
    central_longitude=pdef.slon,
    standard_parallels=(pdef.Struelat, pdef.Ntruelat),
    false_easting=pdef.iref * pdef.dx,
    false_northing=pdef.jref * pdef.dy)
# %%

fig, ax = plt.subplots(1, 3, figsize=(19, 4),
                       subplot_kw={'projection': crs_new})
da.land.plot(ax=ax[0], transform=crs_new)
ax[0].title.set_text('crs_new')
da.land.plot(ax=ax[1], transform=crs)
ax[1].title.set_text('crs')
da = da.set_coords(['XLONG', 'XLAT'])
da.land.plot.pcolormesh(
    ax=ax[2], transform=ccrs.PlateCarree(), x="XLONG", y="XLAT")
ax[0].plot(0, 0, marker="o", markersize=3)
ax[1].plot(0, 0, marker="o", markersize=3)
ax[2].plot(0, 0, marker="o", markersize=3)
ax[2].title.set_text('LatLon')
plt.suptitle('map crs_new')
plt.savefig('test.png',
            dpi=150, bbox_inches='tight')

# %%
fig, ax = plt.subplots(1, 3, figsize=(19, 4),
                       subplot_kw={'projection': crs_new2})
da.land.plot(ax=ax[0], transform=crs_new2)
ax[0].title.set_text('crs_new2')
da.land.plot(ax=ax[1], transform=crs)
ax[1].title.set_text('crs')
da = da.set_coords(['XLONG', 'XLAT'])
da.land.plot.pcolormesh(
    ax=ax[2], transform=ccrs.PlateCarree(), x="XLONG", y="XLAT")
ax[0].plot(0, 0, marker="o", markersize=3)
ax[1].plot(0, 0, marker="o", markersize=3)
ax[2].plot(0, 0, marker="o", markersize=3)
ax[2].title.set_text('LatLon')
plt.suptitle('map crs_new2')
plt.savefig('test2.png',
            dpi=150, bbox_inches='tight')


# %%
'''
import numpy as np
# caculate numbers

lon2d = da.XLONG.values
lat2d = da.XLAT.values

ycoord = np.linspace(0, (pdef.jsize - 1) * pdef.dx, pdef.jsize)
xcoord = np.linspace(0, (pdef.isize - 1) * pdef.dy, pdef.isize)
# case 1
crs = get_data_projection(ctl, globe=globe)
wgs_proj = ccrs.Geodetic(globe=globe)
xx, yy = np.meshgrid(xcoord, ycoord)
x, y, _ = crs.transform_points(wgs_proj, lon2d, lat2d).T
x - xx.T

# %%
# case2
wgs_proj = ccrs.Geodetic(globe=globe)
xx, yy = np.meshgrid(xcoord, ycoord)
x, y, _ = crs_new.transform_points(wgs_proj, lon2d, lat2d).T
x - xx.T
'''