import os
import pandas as pd
import numpy as np
from logging import debug
import matplotlib
# matplotlib.use('Agg')
import matplotlib.pyplot as plt
from ....transforms.basic.regrid.utils.reproject import reproject_emissions
from .....utils.check.errclass import CifError
[docs]
def make_landuse(domain):
dir_landuse = domain.dir_landuse
if not os.path.isdir(dir_landuse):
raise CifError(
f"Could not initialize landuse as dir_landuse = '{dir_landuse}' is not a valid path. Please check your yaml")
# Target domain
zlon = domain.zlon
zlat = domain.zlat
zlonc = domain.zlonc
zlatc = domain.zlatc
nlon = domain.nlon
nlat = domain.nlat
# Define GLCF domain
nlatbig = 21600
nlonbig = 43200
latmax = 89.995833333333333
lonmin = -179.9958333333333333
zinter = 0.008333333333333333
latbig = latmax - zinter * np.arange(nlatbig)
lonbig = lonmin + zinter * np.arange(nlonbig)
# Check what GLCF window is needed
ilatmax_all = np.argmin(np.abs(latbig - zlatc.max()))
ilatmin_all = np.argmin(np.abs(latbig - zlatc.min()))
ilonmin_all = np.argmin(np.abs(lonbig - zlonc.min()))
ilonmax_all = np.argmin(np.abs(lonbig - zlonc.max()))
# Initialize GLCF array
debug(f"Initialize empty GLCF array of size ({ilatmin_all - ilatmax_all}x{ilonmax_all - ilonmin_all})")
glcf_data = np.zeros((ilatmin_all - ilatmax_all, ilonmax_all - ilonmin_all))
nlat_out, nlon_out = glcf_data.shape
# Extract the data
debug("Reading GLCF data")
zones = ["ea", "na", "sa", "af", "ap"]
for zone in zones:
debug(f"Reading zone: {zone}")
file_zone = f"{dir_landuse}/{zone}0500ag.asc"
if not os.path.isfile(file_zone):
raise CifError(f"Could not extract GLCF from {file_zone}. The file is not here.")
with open(file_zone, "r") as f:
ncols = int(f.readline().split()[1])
nrows = int(f.readline().split()[1])
xllcorner = float(f.readline().split()[1])
yllcorner = float(f.readline().split()[1])
cellsize = float(f.readline().split()[1])
nodata = f.readline()
lonzone = xllcorner + cellsize * np.arange(ncols)
latzone = yllcorner + cellsize * np.arange(nrows - 1, -1, -1)
# Where is the domain in the zone
ilatmax_dom = max(0, int(np.round((zlatc.max() - latzone.max()) / cellsize)))
ilatmin_dom = min(nlat_out, int(np.round((zlatc.max() - latzone.min()) / cellsize)))
ilonmin_dom = max(0, int(np.round((lonzone.min() - zlonc.min()) / cellsize)))
ilonmax_dom = min(nlon_out, int(np.round((lonzone.max() - zlonc.min()) / cellsize)))
ilatmax_dom = int(np.round((zlatc.max() - latzone.max()) / cellsize))
ilatmin_dom = int(np.round((zlatc.max() - latzone.min()) / cellsize))
ilonmin_dom = int(np.round((lonzone.min() - zlonc.min()) / cellsize))
ilonmax_dom = int(np.round((lonzone.max() - zlonc.min()) / cellsize))
print(ilatmax_dom, ilatmin_dom, ilonmax_dom, ilonmin_dom)
print(latzone.min(), latzone.max(), lonzone.min(), lonzone.max())
# Stop here if domain is not in zone
# if ilatmin_dom - ilatmax_dom <= 1 or ilonmax_dom - ilonmin_dom <= 1:
if ilatmax_dom >= nlat_out or ilatmin_dom <= 0 \
or ilonmax_dom <= 0 or ilonmin_dom >= nlon_out:
debug(f"Target domain is not in zone {zone}")
continue
# Read the data
for ilat in range(nrows):
ln = f.readline()
if 0 <= ilat + ilatmax_dom < nlat_out:
print(ilat)
fland = np.array(
list(map(float, ln.split()))
)[max(0, -ilonmin_dom):
max(0, -ilonmin_dom)
+ min(nlon_out, ilonmax_dom)
- max(0, ilonmin_dom)
]
mask = glcf_data[
ilat + ilatmax_dom,
max(0, ilonmin_dom):min(nlon_out, ilonmax_dom)
] == 0
glcf_data[ilat + ilatmax_dom,
max(0, ilonmin_dom):min(nlon_out, ilonmax_dom)
][mask] += fland[mask]
# Plot intermediate map in case
plt.figure(figsize=(21, 11))
plt.imshow(glcf_data)
plt.savefig(f"{domain.workdir}/domain/LANDUSE/map_check.png")
plt.close()
# Now reproject GLCF to domain
zlonc_in = np.concatenate([lonbig - zinter / 2, [lonbig[-1] + zinter / 2]])
zlatc_in = np.concatenate([latbig + zinter / 2, [latbig[-1] - zinter / 2]])
zlatc_in, zlonc_in = np.meshgrid(zlatc_in[ilatmax_all:ilatmin_all + 1],
zlonc_in[ilonmin_all:ilonmax_all + 1],
indexing="ij")
zlatc_in = np.flip(zlatc_in, axis=0)
glcf_data = np.flip(glcf_data, axis=0)
nlandin = 14
weights = reproject_emissions(
0 * zlonc_in[1:, 1:],
zlonc_in,
zlatc_in,
zlonc,
zlatc,
orig_regular=True,
orig_unstructured=False,
return_weight=True
)
flandin = np.array([glcf_data == ilan for ilan in range(nlandin)])
fland = np.zeros((nlandin, nlat, nlon))
iout, jout = np.unravel_index(range(nlat * nlon),
(nlat, nlon), order="F")
for iland in range(nlandin):
fland[iland, iout, jout] = np.nansum(
flandin[iland, weights["i"], weights["j"]] * weights["wgt"], axis=-1)
# Now aggregate landuse
nlandout = 9
file_aggregation = domain.file_aggregation
if not os.path.isfile(file_aggregation):
raise CifError(f"Could not find file_aggregation (LANDUSE_AGGREGATION): {file_aggregation}")
aggreg_prop = pd.read_csv(file_aggregation, sep=r"\s+", usecols=range(9))
flandout = np.sum(
fland[:, np.newaxis, :, :] * aggreg_prop.values[:, :, np.newaxis, np.newaxis],
axis=0)
# Visualize output before saving
for ilan in range(nlandout):
plt.figure(figsize=(21, 11))
plt.imshow(np.flip(flandout[ilan], axis=0))
plt.savefig(f"{domain.workdir}/domain/LANDUSE/map_out_LU_{ilan}.png")
plt.close()
# Dump to final file
flu = f"{domain.workdir}/domain/LANDUSE/LANDUSE_{domain.domid}"
np.savetxt(flu, flandout.reshape((nlandout, -1)).T,
fmt='%.3f', delimiter=" ")
