"""
Reproject data
==============

In this example we will see how to reproject vector and raster datasets.

"""
import matplotlib.pyplot as plt
import pyproj

# %%
# We'll start with the imports
import imod

# sphinx_gallery_thumbnail_number = -1

# %%
# Reproject points
# ----------------

wgs84 = pyproj.CRS("EPSG:4326")
amersfoort = pyproj.CRS("EPSG:28992")
lon = [5.053, 4.479, 5.722]
lat = [52.201, 52.009, 52.19246]
x, y = pyproj.transform(wgs84, amersfoort, lat, lon)

print(x, y)

# %%
# Reproject vector dataset
# ------------------------
#
# In this case, the shapefile is imported using geopandas, obtaining a
# GeoDataFrame. GeoPandas has the option geopandas.GeoSeries.to_crs to directly
# change the coordinate system of a geopandas GeoDataFrame.

temp_dir = imod.util.temporary_directory()

lakes = imod.data.lakes_shp(temp_dir)
print(lakes.crs)

lakes_wgs84 = lakes.to_crs(epsg=4326)
print(lakes_wgs84.crs)

# %%
# Reproject raster dataset
# ------------------------
#
# imod-python has the function :py:func:`imod.prepare.reproject`. There are
# three options:
#
# 1. Resample to a new cellsize
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Use the same projection: provide only "like".
#
# Importing the original file, which has a cellsize of 100.0 m and its EPSG is 28992:

ahn = imod.data.ahn()["ahn"]

fig, ax = plt.subplots()
ahn.plot.imshow(ax=ax)

# %%
# We'll create our like grid Resampling DataArray to a new cellsize of 50.0 m,
# by creating a like DataArray first:
xmin = 90950
xmax = 115650
ymax = 467550
ymin = 445850
cellsize = 50.0
dx = cellsize
dy = -cellsize

like = imod.util.empty_2d(dx, xmin, xmax, dy, ymin, ymax)

# %%
# Apply the :py:func:`imod.prepare.reproject` function. The new dataset will
# have a 50 m resolution:
ahn_50m = imod.prepare.reproject(source=ahn, like=like, method="average")
print(ahn_50m.res)

# %%
# 2. Only reproject
# ~~~~~~~~~~~~~~~~~
#
# Only provide the source coordinate reference system (src_crs) and the target
# coordinate reference system (dst_crs). In this case, to reproject from
# EPSG:28992 to EPSG:32631:

ahn_utm = imod.prepare.reproject(source=ahn, src_crs="EPSG:28992", dst_crs="EPSG:32631")
print(ahn_utm.res)

# %%
# If we plot, notice that the grid is slightly "rotated". This is caused by the
# reprojection.
fig, ax = plt.subplots()
ahn_utm.plot.imshow(ax=ax)

# %%
# 3. Reproject and resample to a specific domain
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# Provide "src_crs", "dst_crs" and "like".
# The resulting dataset will have a cellsize of 50m and
# it's coordinate system will be EPSG:32631:

ahn_utm_50m = imod.prepare.reproject(
    source=ahn, like=like, src_crs="EPSG:28992", dst_crs="EPSG:32631"
)
print(ahn_utm_50m.res)
print(ahn_utm_50m.crs)