from osgeo import gdal, ogr, osr
import pandas as pd
from logging import info, debug
import numpy as np
import itertools
from ......utils.check.errclass import CifError, CifValueError
new_pyproj = True
try:
from pyproj import Proj, Transformer
except ImportError:
from pyproj import Proj, transform
new_pyproj = False
[docs]
def reproject_emissions(
emis_orig,
zlonc_orig,
zlatc_orig,
zlonc_target,
zlatc_target,
resol=10,
option="mean",
wk_proj="wgs84",
orig_regular=True,
return_weight=False,
orig_unstructured=False,
orig_lon_cyclic=False,
target_unstructured=False,
rounding_domain=8
):
# Calculate the number of bands
if len(emis_orig.shape) == 2:
debug("Reprojecting 2D data")
nbands = 1
emis_orig = emis_orig[np.newaxis, ...]
elif len(emis_orig.shape) == 3:
debug("Reprojecting 3D data")
nbands = emis_orig.shape[0]
debug(f"Number of bands: {nbands}")
# Shifting longitudes beyond 180 degree east
# Unwrap and makes sure that target is in same modulo than orig
if wk_proj == "wgs84":
zlonc_orig = (zlonc_orig + 180) % 360 - 180
zlonc_target = (zlonc_target + 180) % 360 - 180
# Unwrap
if not orig_unstructured:
zlonc_orig = np.degrees(np.unwrap(np.radians(zlonc_orig)))
else:
zlonc_orig = np.degrees(np.unwrap(np.radians(zlonc_orig), axis=0))
if not target_unstructured:
zlonc_target = np.degrees(np.unwrap(np.radians(zlonc_target)))
else:
zlonc_target = np.degrees(
np.unwrap(np.radians(zlonc_target), axis=0))
# Put to the same range of degrees
zlonc_target += 360 * (
np.median(np.floor((zlonc_orig - 180) / 360))
- np.median(np.floor((zlonc_target - 180) / 360))
)
# Projection definition
wgs84_pyproj = "epsg:4326"
if wk_proj == "wgs84":
ref_proj = osr.SpatialReference()
ref_proj.ImportFromProj4(
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
)
ref_pyproj = "epsg:4326"
else:
ref_proj = osr.SpatialReference()
ref_proj.ImportFromProj4(wk_proj)
ref_pyproj = wk_proj
# Convert lon/lat domain coordinates to reference coordinate system
if new_pyproj:
t = Transformer.from_crs(wgs84_pyproj, ref_pyproj, always_xy=True)
xc_orig, yc_orig = t.transform(zlonc_orig, zlatc_orig)
xc_target, yc_target = t.transform(zlonc_target, zlatc_target)
else:
xc_orig, yc_orig = transform(
Proj(init=wgs84_pyproj), Proj(init=ref_pyproj),
zlonc_orig, zlatc_orig)
xc_target, yc_target = transform(
Proj(init=wgs84_pyproj), Proj(init=ref_pyproj),
zlonc_target, zlatc_target)
# Domains dimensions
nmerid_orig, nzonal_orig = zlonc_orig[:-1, :-1].shape
nmerid_target, nzonal_target = zlonc_target[:-1, :-1].shape
# Check that the dimensions of emissions are compatible with coordinates
if not orig_unstructured:
nmerid_emis, nzonal_emis = emis_orig.shape[1:]
if nmerid_emis != nmerid_orig or nzonal_emis != nzonal_orig:
raise CifError(
"Warning: Emission data shape is "
"not compatible of that of coordinates. \n"
"Consider transposing, or shake how the data is read. \n"
f"Emission shape: {emis_orig.shape}\n"
f"Coordinate shape: {zlonc_orig[:-1, :-1].shape}")
else:
nmerid_emis, nzonal_emis = emis_orig.shape[1:]
if nmerid_emis != 1 and nmerid_emis != zlonc_orig[0].size:
raise CifError(
"Warning: Emission data shape is "
"not compatible of that of coordinates. "
"Consider transposing, or shake how the data is read. . \n"
f"Emission shape: {emis_orig.shape}\n"
f"Coordinate shape: {zlonc_orig[0].size}")
# TODO: make it work with unstructured domains
raise CifError(
"Unstructured domains in origin data not yet implemented in reproject")
# Check that the original domain is really regular
if orig_regular:
if len(np.unique(np.diff(xc_orig, axis=0).round(3))) > 1:
orig_regular = False
if len(np.unique(np.diff(xc_orig, axis=1).round(3))) > 1:
orig_regular = False
if len(np.unique(np.diff(yc_orig, axis=0).round(3))) > 1:
orig_regular = False
if len(np.unique(np.diff(yc_orig, axis=1).round(3))) > 1:
orig_regular = False
# Rasterize to regular grid the origin emissions if necessary
xpxl_ref = None
if not orig_regular:
debug("Irregular domain: Reprojecting original domain "
"to a finer regular domain")
# Rasterize to regular grid the origin emissions if necessary
debug("Vectorizing the original domain first")
(
vector_grid_driver,
vector_grid_ds,
vector_grid_layer,
) = domain2polygons(xc_orig, yc_orig, wk_proj, data=emis_orig)
# Rasterize irregular origin data
x_min, x_max, y_min, y_max = vector_grid_layer.GetExtent()
if not target_unstructured:
pixel_size = (
max(
np.median(np.ediff1d(xc_target)),
np.median(np.ediff1d(xc_target.T)),
)
/ resol
)
else:
pixel_size = (
(xc_target.max() - xc_target.min())
/ np.sqrt(xc_target.shape[1])
) / resol
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
geotransform = (x_min, pixel_size, 0, y_max, 0, -pixel_size)
output_raster = gdal.GetDriverByName("MEM").Create(
"", x_res, y_res, nbands + 2, gdal.GDT_Float32
)
output_raster.SetGeoTransform(geotransform)
output_raster.SetProjection(ref_proj.ExportToWkt())
debug("Rasterize vector domain to finer scale:")
debug("Target raster charecteristics:")
debug(f"xmin/xmax/ymin/ymax: {x_min}/{x_max}/{y_min}/{y_max}")
debug(f"resolution x/y: {pixel_size}")
debug(f"Raster x/y size: {x_res}/{y_res}")
for iband in range(nbands):
gdal.RasterizeLayer(
output_raster,
[iband + 1],
vector_grid_layer,
None,
None,
[0],
options=[f"ATTRIBUTE=field_{iband}"],
)
# Rasterize i/j
gdal.RasterizeLayer(
output_raster,
[nbands + 1],
vector_grid_layer,
None,
None,
[0],
options=["ATTRIBUTE=index_i"],
)
gdal.RasterizeLayer(
output_raster,
[nbands + 2],
vector_grid_layer,
None,
None,
[0],
options=["ATTRIBUTE=index_j"],
)
else:
dx = np.unique(np.round(np.diff(xc_orig, axis=1), 6))
if wk_proj == "wgs84":
# Unwrap geometry
xc_bis = np.degrees(np.unwrap(np.radians(xc_orig)))
dx = np.unique(np.round(np.diff(xc_bis, axis=1), rounding_domain))
xc_bis -= 360 * np.round(xc_bis.mean() / 360)
xc_orig = xc_bis[:]
if np.size(dx) > 1:
raise CifValueError(
f"The domain is not regular.\n"
f"The list of possible deltas is: {dx}.\n"
"If all delta are very similar apart from rounding difference, "
"it is possible to force a constant delta by using the parameter "
f"'rounding_domain' (={rounding_domain} presently) in the yml. "
"Please check the documentation for "
"the transform plugin 'regrid/std'"
)
raise CifError(
f"Warning! The domain is not regular. I found the following possible dx values: {dx}")
dx = dx[0]
# Shift original array spanning beyond 180 degrees east
if wk_proj == "wgs84":
if xc_orig[0, 0] + (nzonal_orig - 1) * dx > 180:
xpxl_ref = int(
np.round(
(180 - xc_orig[0, 0]) / dx
)
)
xc_orig = np.concatenate(
(xc_orig[:, xpxl_ref:], xc_orig[:, :xpxl_ref]), axis=1
)
yc_orig = np.concatenate(
(yc_orig[:, xpxl_ref:], yc_orig[:, :xpxl_ref]), axis=1
)
emis_orig = np.concatenate(
(emis_orig[..., xpxl_ref:], emis_orig[..., :xpxl_ref]),
axis=2
)
# Deal with poles
geotransform = (
xc_orig[0, 0],
dx,
0,
yc_orig[0, 0],
0,
yc_orig[1, 1] - yc_orig[0, 0],
)
# Create Raster with all emissions bands
output_raster = gdal.GetDriverByName("MEM").Create(
"", nzonal_orig, nmerid_orig, nbands, gdal.GDT_Float32
)
output_raster.SetGeoTransform(geotransform)
output_raster.SetProjection(ref_proj.ExportToWkt())
# Loop over month to fill the raster
for iband in range(nbands):
# Writes my array to the raster
output_raster.GetRasterBand(iband + 1).WriteArray(emis_orig[iband])
# Create polygons from target domain
vector_grid_driver, vector_grid_ds, vector_grid_layer = domain2polygons(
zlonc_target, zlatc_target, wk_proj, is_regular=not target_unstructured
)
# Compute projection
emis_target = loop_zonal_stats(
vector_grid_layer,
output_raster,
resol=resol,
option=option,
return_weight=return_weight,
orig_lon_cyclic=orig_lon_cyclic
)
# Return only weights if required
if return_weight:
if xpxl_ref is not None:
for k, wgt in enumerate(emis_target):
if wgt[1] == []:
continue
mask = wgt[1] >= xpxl_ref
emis_target[k][1][mask] -= xpxl_ref
emis_target[k][1][~mask] += xpxl_ref
# Reshape weights as dictionary
weights = {
"i": [e[0] for e in emis_target],
"j": [e[1] for e in emis_target],
"wgt": [e[2] for e in emis_target]
}
# Deal with non regular original domain
if not orig_regular:
debug("Reshaping refined original domain to real one when irregular")
bandi = output_raster.GetRasterBand(nbands + 1).ReadAsArray() - 100
bandj = output_raster.GetRasterBand(nbands + 2).ReadAsArray() - 100
out_weights = {"i": [], "j": [], "wgt": []}
ind_target = 0
nmerid_loop = nmerid_target if not target_unstructured else 1
for ii, jj, ww in zip(weights["i"], weights["j"], weights["wgt"]):
j = ind_target // nmerid_loop
ind_target += 1
if j % int(nzonal_target / 10) == 0:
debug(f"Grid cell: {j} out of {nzonal_target}")
tmp_i = bandi[ii, jj][:, np.newaxis]
tmp_j = bandj[ii, jj][:, np.newaxis]
unique, indices_target = np.unique(
np.concatenate([tmp_i, tmp_j], axis=1),
axis=0, return_inverse=True)
wgt_target = np.zeros(len(unique))
np.add.at(wgt_target, indices_target, ww)
out_weights["i"].append(unique[:, 0].flatten().astype(int))
out_weights["j"].append(unique[:, 1].flatten().astype(int))
out_weights["wgt"].append(wgt_target)
weights = out_weights
# Converts weights and i/j to a matrix padded with nans for varying
# lengths of i/j/weights
debug("Reshaping outputs into a full shape matrix before return weights")
weights = {
"i": np.array(list(
itertools.zip_longest(*weights["i"], fillvalue=0))).T,
"j": np.array(list(
itertools.zip_longest(*weights["j"], fillvalue=0))).T,
"wgt": np.array(list(
itertools.zip_longest(*weights["wgt"], fillvalue=np.nan))).T
}
return weights
emis_target = np.reshape(emis_target, (nmerid_target, nzonal_target, -1))
emis_target = np.transpose(emis_target, axes=(2, 1, 0))
# Removes NaNs
emis_target[np.isnan(emis_target)] = 0.0
# if only one band, returns only the 2D dataset,
# otherwise, returns everything
if nbands == 1:
return emis_target[0]
else:
return emis_target
[docs]
def domain2polygons(zlonc_target, zlatc_target,
wk_proj="wgs84", data=None, is_regular=True):
debug(
f"Vectorizing the domain of shape {zlonc_target.shape} and is_regular={is_regular}")
debug("This can take a while...")
if is_regular:
nmerid_target, nzonal_target = zlonc_target[:-1, :-1].shape
else:
nmerid_target = 1
nzonal_target = zlonc_target.shape[1]
if data is not None:
nbands = data.shape[0]
# Reference projections
wgs84_pyproj = "epsg:4326"
if wk_proj == "wgs84":
ref_proj = osr.SpatialReference()
ref_proj.ImportFromProj4(
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
)
ref_pyproj = "epsg:4326"
else:
ref_proj = osr.SpatialReference()
ref_proj.ImportFromProj4(wk_proj)
ref_pyproj = wk_proj
# Convert lon/lat domain coordinates to reference coordinate system
if new_pyproj:
t = Transformer.from_crs(wgs84_pyproj, ref_pyproj, always_xy=True)
xc, yc = t.transform(zlonc_target, zlatc_target)
else:
xc, yc = transform(
Proj(init=wgs84_pyproj), Proj(init=ref_pyproj),
zlonc_target, zlatc_target)
# define the layer
vector_grid_driver = ogr.GetDriverByName("MEMORY")
vector_grid_ds = vector_grid_driver.CreateDataSource("memData")
vector_grid_layer = vector_grid_ds.CreateLayer(
"grid", ref_proj, geom_type=ogr.wkbPolygon
)
idField = ogr.FieldDefn("index_i", ogr.OFTInteger)
vector_grid_layer.CreateField(idField)
idField = ogr.FieldDefn("index_j", ogr.OFTInteger)
vector_grid_layer.CreateField(idField)
if data is not None:
for iband in range(nbands):
idField = ogr.FieldDefn(f"field_{iband}", ogr.OFTReal)
vector_grid_layer.CreateField(idField)
featureDefn = vector_grid_layer.GetLayerDefn()
for j in range(nzonal_target):
if j % int(nzonal_target / 10) == 0:
debug(f"Vectorizing: {j} out of {nzonal_target}")
for i in range(nmerid_target):
if is_regular:
lon_tmp = [
xc[i, j],
xc[i, j + 1],
xc[i + 1, j + 1],
xc[i + 1, j],
xc[i, j],
]
if wk_proj == "wgs84":
lon_tmp = np.degrees(np.unwrap(np.radians(lon_tmp)))
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
_ = ring.AddPoint(lon_tmp[0], yc[i, j])
_ = ring.AddPoint(lon_tmp[1], yc[i, j + 1])
_ = ring.AddPoint(lon_tmp[2], yc[i + 1, j + 1])
_ = ring.AddPoint(lon_tmp[3], yc[i + 1, j])
_ = ring.AddPoint(lon_tmp[4], yc[i, j])
else:
lon_tmp = xc[:, j]
npoints = len(lon_tmp)
if wk_proj == "wgs84":
lon_tmp = np.degrees(np.unwrap(np.radians(lon_tmp)))
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
for i in range(npoints + 1):
_ = ring.AddPoint(lon_tmp[i % npoints], yc[i % npoints, j])
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
# Add feature to layer
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
# Keep in memory original i/j
feature.SetField("index_i", i + 100)
feature.SetField("index_j", j + 100)
# Add data if any
if data is not None:
for iband in range(nbands):
feature.SetField(
f"field_{iband}", np.double(data[iband, i, j])
)
vector_grid_layer.CreateFeature(feature)
feature = None
debug("Vectorizing is done.")
return vector_grid_driver, vector_grid_ds, vector_grid_layer
####################
# ZONAL STATISTICS #
####################
[docs]
def zonal_stats(feat, raster, resol=10, option="mean", return_weight=False,
orig_lon_cyclic=False, **kwargs):
# Create for target raster the same projection as for the value raster
raster_srs = osr.SpatialReference()
raster_srs.ImportFromWkt(raster.GetProjectionRef())
# Create temporary layer from feature
drv = ogr.GetDriverByName("MEMORY")
ds = drv.CreateDataSource("memData")
layer = ds.CreateLayer("grid", raster_srs, geom_type=ogr.wkbPolygon)
layer.CreateFeature(feat)
# Get raster georeference info
geotransform = raster.GetGeoTransform()
xOrigin = geotransform[0]
yOrigin = geotransform[3]
pixelWidth = geotransform[1]
pixelHeight = geotransform[5]
xWest = min(xOrigin, xOrigin + pixelWidth * raster.RasterXSize)
xEast = max(xOrigin, xOrigin + pixelWidth * raster.RasterXSize)
ySouth = min(yOrigin, yOrigin + pixelHeight * raster.RasterYSize)
yNorth = max(yOrigin, yOrigin + pixelHeight * raster.RasterYSize)
# Get extent of feat
geom = feat.GetGeometryRef()
if geom.GetGeometryName() == "MULTIPOLYGON":
count = 0
pointsX = []
pointsY = []
for _ in geom:
geomInner = geom.GetGeometryRef(count)
ring = geomInner.GetGeometryRef(0)
numpoints = ring.GetPointCount()
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)
count += 1
elif geom.GetGeometryName() == "POLYGON":
ring = geom.GetGeometryRef(0)
numpoints = ring.GetPointCount()
pointsX = []
pointsY = []
for p in range(numpoints):
lon, lat, z = ring.GetPoint(p)
pointsX.append(lon)
pointsY.append(lat)
else:
raise CifError(
"ERROR: Geometry needs to be either Polygon or Multipolygon"
)
xmin = min(pointsX)
xmax = max(pointsX)
ymin = min(pointsY)
ymax = max(pointsY)
if xmax <= xWest or xmin >= xEast or ymax <= ySouth or ymin >= yNorth:
debug("Outside input domain")
debug(geotransform)
debug([raster.RasterXSize, raster.RasterYSize])
debug([xmin, xmax, ymin, ymax])
debug([xWest, xEast, ySouth, yNorth])
debug(pointsX)
debug(pointsY)
if return_weight:
return [], [], []
else:
return raster.RasterCount * [np.nan]
# Specify offset and rows and columns to read
xoff1 = int(np.floor((xmin - xOrigin) / pixelWidth))
xoff2 = int(np.floor((xmax - xOrigin) / pixelWidth))
xoff = max(0, min(xoff1, xoff2))
yoff1 = int(np.floor((ymin - yOrigin) / pixelHeight))
yoff2 = int(np.floor((ymax - yOrigin) / pixelHeight))
yoff = max(0, min(yoff1, yoff2))
xcount = min(max(xoff1, xoff2) + 1, raster.RasterXSize - 1) - xoff + 1
ycount = min(max(yoff1, yoff2) + 1, raster.RasterYSize - 1) - yoff + 1
# Create memory target raster
target_ds = gdal.GetDriverByName("MEM").Create(
"", xcount * resol, ycount * resol, 1, gdal.GDT_Byte
)
target_ds.SetGeoTransform(
(
xOrigin + xoff * pixelWidth,
pixelWidth / resol,
0,
yOrigin + yoff * pixelHeight,
0,
pixelHeight / resol,
)
)
# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[1])
bandmask = target_ds.GetRasterBand(1)
datamask = bandmask.ReadAsArray().astype(float)
datamask2 = [
[datamask[i::resol, j::resol] for j in range(resol)]
for i in range(resol)
]
datamask2 = np.array(datamask2).mean(axis=(0, 1))
# If asks weights only
if return_weight:
meshj, meshi = np.meshgrid(
np.arange(xoff, xoff + xcount), np.arange(yoff, yoff + ycount)
)
mask = datamask2 > 0
return meshi[mask], meshj[mask], datamask2[mask] / datamask2.sum()
# Equivalent clipped raster
nbands = raster.RasterCount
raster_locmean = []
for iband in range(1, nbands + 1):
banddataraster = raster.GetRasterBand(iband)
dataraster = banddataraster.ReadAsArray(
xoff, yoff, xcount, ycount
).astype(float)
# Averaging over the area corresponding to the pixel
zoneraster = dataraster * datamask2 / datamask2.sum()
if option == "mean":
raster_locmean.append(np.sum(zoneraster))
elif option == "median":
raster_locmean.append(np.median(zoneraster))
elif option == "min":
raster_locmean.append(zoneraster[np.argmin(datamask2)])
elif option == "max":
raster_locmean.append(
dataraster[
np.unravel_index(datamask2.argmax(), datamask2.shape)
]
)
else:
raster_locmean.append(np.nan)
# Calculate statistics of zonal raster
return raster_locmean
[docs]
def loop_zonal_stats(
lyr, raster, resol=10, option="mean", return_weight=False,
orig_lon_cyclic=False
):
featList = range(lyr.GetFeatureCount())
statList = []
debug("Looping on target grid cells and finding proportion from "
"original domain")
for FID in featList:
try:
feat = lyr.GetFeature(FID)
meanValue = zonal_stats(
feat,
raster,
resol=resol,
option=option,
return_weight=return_weight,
FID=FID,
orig_lon_cyclic=orig_lon_cyclic
)
except Exception as e:
info(e)
meanValue = raster.RasterCount * [np.nan]
raise e
if FID % int(len(featList) / 10) == 0:
debug(f"Target grid cell: {FID} out of {len(featList)}")
statList.append(meanValue)
if return_weight:
return statList
else:
return np.array(statList)
