Source code for wolfhece.wolf_array._header_wolf
"""
Header class for WOLF arrays.
Contains:
- WOLF_ARRAY_* type constants
- getkeyblock / decodekeyblock helper functions
- header_wolf class (spatial header: origin, resolution, shape, translation, multi-block support)
This module is intentionally free of wx, OpenGL, and heavy GUI dependencies so that
it can be imported cheaply in non-interactive contexts (CLI tools, unit tests, etc.).
Author: HECE - University of Liege, Pierre Archambeau
Date: 2024
Copyright (c) 2024 University of Liege. All rights reserved.
This script and its content are protected by copyright law. Unauthorized
copying or distribution of this file, via any medium, is strictly prohibited.
"""
import os
import json
import logging
import numpy as np
from pathlib import Path
from typing import Union, Literal
from copy import deepcopy
try:
from osgeo import gdal, osr
gdal.UseExceptions()
except ImportError as e:
print(e)
raise Exception('Error importing GDAL library')
try:
from ..PyTranslate import _
except ImportError:
from shapely.geometry import Polygon
# Lazy imports for type hints only
from typing import TYPE_CHECKING
if TYPE_CHECKING:
from ..PyVertexvectors import vector, zone
# ======================================================================
# WOLF array type constants
# ======================================================================
# These types are not supported in Fortran, but they are supported in Python.
# ======================================================================
# Block key helpers
# ======================================================================
[docs]
def getkeyblock(i, addone=True) -> str:
"""
Name/Key of a block in the dictionnary of a WolfArrayMB instance
For Fortran compatibility, addone is True by default so first block is "block1" and not "block0"
"""
if addone:
return 'block' + str(i + 1)
else:
return 'block' + str(i)
[docs]
def decodekeyblock(key, addone=True) -> int:
"""
Decode key of a block in the dictionnary of a WolfArrayMB instance
For Fortran compatibility, addone is True by default so first block is "block1" and not "block0"
"""
if addone:
return int(key[5:])
else:
return int(key[5:]) - 1
# ======================================================================
# header_wolf
# ======================================================================
[docs]
class header_wolf():
"""
Header of WolfArray
In case of a mutliblock, the header have informations about all the blocks in head_blocks dictionnary.
Block keys are generated by "getkeyblock" function
"""
# FIXME It'd be wise to put the multiblock case into another class.
# for example "header_wolf_MB" else one could construct hierearchies
# of headers which don't exist in practice.
def __init__(self) -> None:
"""
Origin (origx, origy, [origz]) is the point in local space from which every other coordinates are measured.
Translation (translx, transly, [translz]) is the translation of the origin in global space. If translation is null, the origin is the same in local and global space. :-)
Resolution (dx, dy, [dz]) is the spatial resolution of the array.
Nullvalue is the value of the null value in the array.
(nbx, nby, [nbz]) are the number of cells in the array along X and Y [and Z]. It is the shape of the array.
@property nbdims is the number of dimensions of the array (2 or 3)
"""
self.head_blocks = {}
def __str__(self) -> str:
""" Return a string representation of the header """
ret = ''
ret += _('Shape : {} x {} \n').format(self.nbx, self.nby)
ret += _('Resolution : {} x {} \n').format(self.dx, self.dy)
ret += _('Spatial extent : \n')
ret += _(' - Origin : ({} ; {}) \n').format(self.origx, self.origy)
ret += _(' - End : ({} ; {}) \n').format(self.origx + self.nbx * self.dx, self.origy +self.nby * self.dy)
ret += _(' - Width x Height : {} x {} \n').format(self.nbx * self.dx, self.nby * self.dy)
ret += _(' - Translation : ({} ; {})\n').format(self.translx, self.transly)
ret += _('Null value : {}\n\n'.format(self.nullvalue))
if len(self.head_blocks) > 0:
ret += _('Number of blocks : {}\n\n').format(len(self.head_blocks))
for key, value in self.head_blocks.items():
ret += _('Block {} : \n\n').format(key)
ret += str(value)
return ret
@property
@epsg.setter
def epsg(self, value:int | str):
if isinstance(value, str):
if value in ['Belgique 1972', 'EPSG:31370', '31370', 'Belgium 1972']:
self._epsg = 31370
elif value in ['RGF93 / Lambert-93', 'EPSG:2154', '2154', 'France Lambert-93']:
self._epsg = 2154
elif value in ['Detutschland GK3', 'EPSG:31467', '31467', 'Germany GK3', 'Allemagne GK3']:
self._epsg = 31467
else:
try:
epsg_code = int(value.split(':')[1])
self._epsg = epsg_code
except (IndexError, ValueError):
logging.warning(_('Could not parse EPSG code from string: {}').format(value))
self._epsg = None
elif isinstance(value, int):
self._epsg = value
@property
@nullvalue.setter
def nullvalue(self, value:float):
self._nullvalue = value
@property
[docs]
def nbdims(self):
if self.nbz == 0:
if self.nbx > 0 and self.nby > 0:
return 2
else:
return 0
elif self.nbz > 0:
return 3
else:
raise Exception(_('The number of dimensions is not correct'))
@nbdims.setter
def nbdims(self, value):
logging.warning(_('nbdims was an attribute of header_wolf.\nIt is now a read-only property.\nPlease use nbx, nby and nbz instead to define the shape of the array'))
raise Exception(_('This property is read-only'))
@property
[docs]
def shape(self):
if self.nbdims == 2:
return (self.nbx, self.nby)
elif self.nbdims == 3:
return (self.nbx, self.nby, self.nbz)
else:
return (0, 0)
@shape.setter
def shape(self, value:tuple[int]):
if len(value) == 3:
self.nbx = value[0]
self.nby = value[1]
self.nbz = value[2]
elif len(value) == 2:
self.nbx = value[0]
self.nby = value[1]
self.nbz = 0
else:
raise Exception(_('The number of dimensions is not correct'))
@property
def __getitem__(self, key:Union[int,str]=None):
"""
Return block header
:param key: block's index (0-based) or key (str)
:return: header_wolf instance if key is found, None otherwise
"""
if key is None:
return self
if isinstance(key,int):
_key = getkeyblock(key)
else:
_key = key
if _key in self.head_blocks.keys():
return self.head_blocks[_key]
else:
return None
def __setitem__(self, key:Union[int,str], value:"header_wolf"):
"""
Set block header
:param value: tuple (key, header_wolf)
'key' can be an int (0-based) or a str
If str, please use getkeyblock function to generate the key
"""
if isinstance(key,int):
_key = getkeyblock(key)
else:
_key = key
self.head_blocks[_key] = deepcopy(value)
@property
@resolution.setter
def resolution(self, value:tuple[float]):
if len(value) == 2:
self.set_resolution(value[0], value[1])
elif len(value) == 3:
self.set_resolution(value[0], value[1], value[2])
@property
@origin.setter
def origin(self, value:tuple[float]):
if len(value) == 2:
self.set_origin(value[0], value[1])
elif len(value) == 3:
self.set_origin(value[0], value[1], value[2])
@property
@translation.setter
def translation(self, value:tuple[float]):
if len(value) == 2:
self.set_translation(value[0], value[1])
elif len(value) == 3:
self.set_translation(value[0], value[1], value[2])
[docs]
def set_resolution(self, dx:float, dy:float, dz:float= 0.):
"""
Set resolution
"""
self.dx = dx
self.dy = dy
self.dz = dz
[docs]
def set_origin(self, x:float, y:float, z:float = 0.):
"""
Set origin
:param x: origin along X
:param y: origin along Y
:param z: origin along Z
"""
self.origx = x
self.origy = y
self.origz = z
[docs]
def get_translation(self):
"""
Return translation
"""
return (self.translx, self.transly, self.translz)
[docs]
def set_translation(self, tr_x:float, tr_y:float, tr_z:float= 0.):
"""
Set translation
:param tr_x: translation along X
:param tr_y: translation along Y
:param tr_z: translation along Z
"""
self.translx = tr_x
self.transly = tr_y
self.translz = tr_z
[docs]
def get_bounds(self, abs=True):
"""
Return bounds in coordinates
:param abs: if True, add translation to (x, y) (coordinate to global space)
:return: tuple of two lists of two floats - ([xmin, xmax],[ymin, ymax])
"""
if abs:
return ([float(self.origx + self.translx), float(self.origx + self.translx + float(self.nbx) * self.dx)],
[float(self.origy + self.transly), float(self.origy + self.transly + float(self.nby) * self.dy)])
else:
return ([float(self.origx), float(self.origx + float(self.nbx) * self.dx)],
[float(self.origy), float(self.origy + float(self.nby) * self.dy)])
[docs]
def get_extent(self, abs=True):
"""
Return extent in coordinates
:param abs: if True, add translation to (x, y) (coordinate to global space)
:return: list of four floats - [xmin, xmax, ymin, ymax]
"""
bounds = self.get_bounds(abs)
return [bounds[0][0], bounds[0][1], bounds[1][0], bounds[1][1]]
[docs]
def get_scale_yoverx(self):
"""
Return the scale of the array (height over width)
"""
if self.dx == 0 or self.dy == 0 or self.nbx == 0 or self.nby == 0:
return 1.
return float((self.nby * self.dy) / (self.nbx * self.dx))
[docs]
def get_bounds_ij(self, abs=False):
"""
Return bounds in indices
Firstly, get_bounds is called to get bounds in coordinates and then get_ij_from_xy is called to get bounds in indices.
:param abs: if True, add translation to (x, y) (coordinate to global space)
"""
mybounds = self.get_bounds(abs)
return (
[self.get_ij_from_xy(mybounds[0][0], mybounds[1][0], abs=abs), self.get_ij_from_xy(mybounds[0][1], mybounds[0][0], abs=abs)],
[self.get_ij_from_xy(mybounds[0][0], mybounds[1][1], abs=abs), self.get_ij_from_xy(mybounds[0][1], mybounds[1][1], abs=abs)])
[docs]
def get_ij_from_xy(self, x:float, y:float, z:float=0., scale:float=1., aswolf:bool=False, abs:bool=True, forcedims2:bool=False) -> Union[tuple[np.int32,np.int32], tuple[np.int32,np.int32,np.int32]]:
"""
Get indices from coordinates
:param x: X coordinate
:param y: Y coordinate
:param z: Z coordinate (optional)
:param scale: scaling of the spatial resolution (dx,dy,[dz])
:param aswolf: if True, return if one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs: if True, remove translation from (x, y, [z]) (coordinate from global space)
:param forcedims2: if True, force to return only 2 indices even if z is supplied
"""
locx = np.float64(x) - self.origx
locy = np.float64(y) - self.origy
locz = np.float64(z) - self.origz
if abs:
locx = locx - self.translx
locy = locy - self.transly
locz = locz - self.translz
i = np.int32(np.floor(locx / (self.dx * scale)))
j = np.int32(np.floor(locy / (self.dy * scale)))
if aswolf:
i += 1
j += 1
if self.nbdims == 3 and not forcedims2:
k = np.int32(np.floor(locz / (self.dz * scale)))
if aswolf:
k += 1
return i, j, k
elif self.nbdims == 2 or forcedims2:
return i, j
[docs]
def get_ij_from_xy_array(self, xy:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True, forcedims2:bool=False) -> np.ndarray:
"""
Get indices from coordinates
:param xy = numpy array containing (x, y, [z]) coordinates - shape (n, 2) or (n, 3)
:param scale = scaling of the spatial resolution (dx,dy,[dz])
:param aswolf = if True, return if one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs = if True, remove translation from (x, y, [z]) (coordinate from global space)
:param forcedims2 = if True, force to return only 2 indices even if z is supplied
:return: numpy array containing (i, j, [k]) indices - shape (n, 2) or (n, 3)
"""
if isinstance(xy,tuple):
if len(xy) == 2:
if (isinstance(xy[0],np.ndarray)) and (isinstance(xy[1],np.ndarray)):
if len(xy[0]) == len(xy[1]):
locxy = np.vstack((xy[0], xy[1])).T
logging.warning(_('get_ij_from_xy_array - xy is a tuple of 2 arrays, it is converted to a 2D array'))
else:
locxy = np.array(xy)
elif isinstance(xy,list):
locxy = np.array(xy)
else:
locxy = xy.copy()
if forcedims2:
locij = np.zeros((locxy.shape[0],2), dtype=np.int32)
else:
locij = np.zeros(locxy.shape, dtype=np.int32)
if locxy.size == 0:
return locij
locxy[:,0] -= self.origx
locxy[:,1] -= self.origy
if abs:
locxy[:,0] -= self.translx
locxy[:,1] -= self.transly
i = np.int32(locxy[:,0] / (self.dx * scale))
j = np.int32(locxy[:,1] / (self.dy * scale))
if aswolf:
i += 1
j += 1
if self.nbdims == 3 and not forcedims2:
locxy[:,2] -= self.origz
if abs:
locxy[:,2] -= self.translz
k = np.int32(locxy[:,2] / (self.dz * scale))
if aswolf:
k += 1
locij[:,0] = i
locij[:,1] = j
locij[:,2] = k
return locij
elif self.nbdims == 2 or forcedims2:
locij[:,0] = i
locij[:,1] = j
return locij
[docs]
def get_xy_from_ij(self, i:int, j:int, k:int=0, scale:float=1., aswolf:bool=False, abs:bool=True) -> Union[tuple[np.float64,np.float64], tuple[np.float64,np.float64,np.float64]]:
"""
Get coordinates from indices
:param i = index along X coordinate
:param j = index along Y coordinate
:param k = index along Z coordinate (optional)
:param scale = scaling of the spatial resolution (dx,dy,[dz])
:param aswolf = if True, input is one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs = if True, add translation to results (x, y, [z]) (coordinate to global space)
"""
i = np.int32(i)
j = np.int32(j)
if aswolf:
# FIXME Put assertion here.
i += -1
j += -1
if abs:
x = (np.float64(i) + .5) * (self.dx * scale) + self.origx + self.translx
y = (np.float64(j) + .5) * (self.dy * scale) + self.origy + self.transly
else:
x = (np.float64(i) + .5) * (self.dx * scale) + self.origx
y = (np.float64(j) + .5) * (self.dy * scale) + self.origy
if self.nbdims == 3:
k = np.int32(k)
if aswolf:
k += -1
if abs:
z = (np.float64(k) - .5) * (self.dz * scale) + self.origz + self.translz
else:
z = (np.float64(k) - .5) * (self.dz * scale) + self.origz
return x, y, z
elif self.nbdims == 2:
return x, y
else:
raise Exception(_("The number of coordinates is not correct"))
[docs]
def transform(self):
""" Return the affine transformation.
Similar to rasterio.transform.Affine
In WOLF, the convention is :
- origin is at the lower-left corner
- the origin is at the corner of the cell dx, dy, so the center of the cell is at dx/2, dy/2
- X axis is along the rows - i index
- Y axis is along the columns - j index
So, the affine transformation is :
(dx, 0, origx + translx + dx /2, 0, dy, origy + transly + dy/2)
"""
from rasterio.transform import Affine
return Affine(self.dx, 0, self.origx + self.translx + self.dx / 2,
0, self.dy, self.origy + self.transly + self.dy / 2)
[docs]
def _transform_gmrio(self):
""" Return the affine transformation.
!! Inverted ij/ji convention !!
Similar to rasterio.transform.Affine
In WOLF, the convention is :
- origin is at the lower-left corner
- the origin is at the corner of the cell dx, dy, so the center of the cell is at dx/2, dy/2
- X axis is along the rows - i index
- Y axis is along the columns - j index
So, the affine transformation is :
(dx, 0, origx + translx + dx /2, 0, dy, origy + transly + dy/2)
"""
from rasterio.transform import Affine
return Affine(0, self.dy, -(self.origy + self.transly),
self.dx, 0, -(self.origx + self.translx))
[docs]
def get_xy_from_ij_array(self, ij:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
"""
Converts array coordinates (numpy cells) to this array's world coodinates.
:param ij = numpy array containing (i, j, [k]) indices - shape (n, 2) or (n, 3)
:param scale = scaling of the spatial resolution (dx,dy,[dz])
:param aswolf = if True, input is one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs = if True, add translation to results (x, y, [z]) (coordinate to global space)
..warning: 'ij' is not the result of np.where() but if you want to use np.where() you can use the following code:
```
np.vstack((ij[0], ij[1])).T
```
"""
if isinstance(ij,tuple):
if len(ij) == 2:
if (isinstance(ij[0],np.ndarray)) and (isinstance(ij[1],np.ndarray)):
if len(ij[0]) == len(ij[1]):
ij = np.vstack((ij[0], ij[1])).T
logging.warning(_('get_xy_from_ij_array - ij is a tuple of 2 arrays, it is converted to a 2D array'))
else:
ij = np.array(ij)
elif isinstance(ij,list):
if len(ij) == 2:
if (isinstance(ij[0],np.ndarray)) and (isinstance(ij[1],np.ndarray)):
if len(ij[0]) == len(ij[1]):
ij = np.vstack((ij[0], ij[1])).T
logging.warning(_('get_xy_from_ij_array - ij is a list of 2 arrays, it is converted to a 2D array'))
else:
ij = np.array(ij)
if abs:
tr_x = self.translx
tr_y = self.transly
tr_z = self.translz
else:
tr_x = 0.
tr_y = 0.
tr_z = 0.
if aswolf:
decali = -1
decalj = -1
decalk = -1
else:
decali = 0
decalj = 0
decalk = 0
xy = np.zeros(ij.shape)
xy[:,0] = (np.float64( (ij[:,0])+decali) + .5) * (self.dx*scale) + self.origx + tr_x
xy[:,1] = (np.float64( (ij[:,1])+decalj) + .5) * (self.dy*scale) + self.origy + tr_y
if self.nbdims == 3 and ij.shape[1]==3:
xy[:,2] = (np.float64( (ij[:,2])+decalk) + .5) * (self.dz*scale) + self.origz + tr_z
return xy
[docs]
def ij2xy(self, i:int, j:int, k:int=0, scale:float=1., aswolf:bool=False, abs:bool=True) -> Union[tuple[np.float64,np.float64], tuple[np.float64,np.float64,np.float64]]:
""" alias for get_xy_from_ij """
return self.get_xy_from_ij(i, j, k, scale, aswolf, abs)
[docs]
def ij2xy_np(self, ij:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
""" alias for get_xy_from_ij_array
:param ij: numpy array containing (i, j, [k]) indices - like np.argwhere() - shape (n, 2) or (n, 3)
:param scale: scaling of the spatial resolution (dx,dy,[dz])
:param aswolf: if True, input is one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs: if True, add translation to results (x, y, [z]) (coordinate to global space)
..warning: 'ij' is not the result of np.where() but if you want to use np.where() you can use the following code:
```
np.vstack((ij[0], ij[1])).T
```
:return: numpy array containing (x, y, [z]) coordinates - shape (n, 2) or (n, 3)
"""
return self.get_xy_from_ij_array(ij, scale, aswolf, abs)
[docs]
def ij2xy_np_ij_from_npwhere(self, ij:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
""" alias for get_xy_from_ij_array but with ij as np.where result
"""
ij = np.vstack((ij[0], ij[1])).T
return self.get_xy_from_ij_array(ij, scale, aswolf, abs)
[docs]
def xy2ij(self, x:float, y:float, z:float=0., scale:float=1., aswolf:bool=False, abs:bool=True, forcedims2:bool=False) -> Union[tuple[np.int32,np.int32], tuple[np.int32,np.int32,np.int32]]:
""" alias for get_ij_from_xy """
return self.get_ij_from_xy(x, y, z, scale, aswolf, abs, forcedims2)
[docs]
def xy2ij_np(self, xy:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
"""
alias for get_ij_from_xy_array
:param xy: numpy array containing (x, y, [z]) coordinates - shape (n, 2) or (n, 3)
:param scale: scaling of the spatial resolution (dx,dy,[dz])
:param aswolf: if True, return if one-based (as Wolf VB6 or Fortran), otherwise 0-based (default Python standard)
:param abs: if True, remove translation from (x, y, [z]) (coordinate from global space)
:param forcedims2: if True, force to return only 2 indices even if z is supplied
:return : numpy array containing (i, j, [k]) indices - shape (n, 2) or (n, 3)
"""
return self.get_ij_from_xy_array(xy, scale, aswolf, abs)
[docs]
def xyz2ijk_np(self, xyz:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
""" alias for get_xy_from_ij_array """
assert xyz.shape[1] == 3, _('xyz must be a 2D array with 3 columns')
return self.get_xy_from_ij_array(xyz, scale, aswolf, abs)
[docs]
def ijk2xyz_np(self, ijk:np.ndarray, scale:float=1., aswolf:bool=False, abs:bool=True) -> np.ndarray:
""" alias for get_xy_from_ij_array """
assert ijk.shape[1] == 3, _('ijk must be a 2D array with 3 columns')
return self.get_xy_from_ij_array(ijk, scale, aswolf, abs)
[docs]
def find_intersection(self, other:"header_wolf", ij:bool = False) -> Union[tuple[list[float],list[float]], tuple[list[list[float]],list[list[float]]]]:
"""
Find the intersection of two header
:param other: other header
:param ij: if True, return indices instead of coordinates
:return: None or tuple of two lists of two floats - ([xmin, xmax],[ymin, ymax]) or indices in each header (if ij=True) [[imin1, imax1], [jmin1, jmax1]], [[imin2, imax2], [jmin2, jmax2]]
"""
mybounds = self.get_bounds()
otherbounds = other.get_bounds()
if otherbounds[0][0] > mybounds[0][1]:
return None
elif otherbounds[1][0] > mybounds[1][1]:
return None
elif otherbounds[0][1] < mybounds[0][0]:
return None
elif otherbounds[1][1] < mybounds[1][0]:
return None
else:
ox = max(mybounds[0][0], otherbounds[0][0])
oy = max(mybounds[1][0], otherbounds[1][0])
ex = min(mybounds[0][1], otherbounds[0][1])
ey = min(mybounds[1][1], otherbounds[1][1])
if ij:
i1, j1 = self.get_ij_from_xy(ox, oy)
i2, j2 = self.get_ij_from_xy(ex, ey)
i3, j3 = other.get_ij_from_xy(ox, oy)
i4, j4 = other.get_ij_from_xy(ex, ey)
return ([[i1, i2], [j1, j2]],
[[i3, i4], [j3, j4]])
else:
return ([ox, ex], [oy, ey])
[docs]
def find_union(self, other:Union["header_wolf", list["header_wolf"]]) -> tuple[list[float],list[float]]:
"""
Find the union of two header
:return: tuple of two lists of two floats - ([xmin, xmax],[ymin, ymax])
"""
if isinstance(other, list):
for cur in other:
assert isinstance(cur, header_wolf), _('All elements in the list must be header_wolf instances')
[ox,ex], [oy,ey] = self.get_bounds()
for cur in other:
otherbounds = cur.get_bounds()
ox = min(ox, otherbounds[0][0])
oy = min(oy, otherbounds[1][0])
ex = max(ex, otherbounds[0][1])
ey = max(ey, otherbounds[1][1])
else:
mybounds = self.get_bounds()
otherbounds = other.get_bounds()
ox = min(mybounds[0][0], otherbounds[0][0])
oy = min(mybounds[1][0], otherbounds[1][0])
ex = max(mybounds[0][1], otherbounds[0][1])
ey = max(mybounds[1][1], otherbounds[1][1])
return ([ox, ex], [oy, ey])
[docs]
def read_txt_header(self, filename:str):
"""
Read informations from header .txt
:param filename: path and filename of the basefile
If filename is a Path object, it is converted to a string
If filename ends with '.tif', nothing is done because infos are in the .tif file
If filename ends with '.flt', a .hdr file must be present and it will be read
Otherwise, a filename.txt file must be present
"""
if isinstance(filename, Path):
filename = str(filename)
locpath = Path(filename)
if filename.endswith('.tif') or filename.endswith('.tiff') :
# from osgeo import gdal
raster:gdal.Dataset
raster = gdal.Open(filename)
geotr = raster.GetGeoTransform()
self.dx = geotr[1]
self.dy = abs(geotr[5])
if geotr[5] < 0.:
self.origx = geotr[0]
self.origy = geotr[3] + raster.RasterYSize * geotr[5]
else:
self.origx = geotr[0]
self.origy = geotr[3]
self.nbx = raster.RasterXSize
self.nby = raster.RasterYSize
"""
https://docs.qgis.org/3.34/en/docs/user_manual/processing_algs/gdal/rasterconversion.html
0 — Use Input Layer Data Type
1 — Byte (Eight bit unsigned integer (quint8))
2 — Int16 (Sixteen bit signed integer (qint16))
3 — UInt16 (Sixteen bit unsigned integer (quint16))
4 — UInt32 (Thirty two bit unsigned integer (quint32))
5 — Int32 (Thirty two bit signed integer (qint32))
6 — Float32 (Thirty two bit floating point (float))
7 — Float64 (Sixty four bit floating point (double))
8 — CInt16 (Complex Int16)
9 — CInt32 (Complex Int32)
10 — CFloat32 (Complex Float32)
11 — CFloat64 (Complex Float64)
12 — Int8 (Eight bit signed integer (qint8))
"""
dtype = raster.GetRasterBand(1).DataType
if dtype == 1:
self.wolftype = WOLF_ARRAY_FULL_INTEGER8
elif dtype in [2,3]:
self.wolftype = WOLF_ARRAY_FULL_INTEGER16
elif dtype in [4,5] :
self.wolftype = WOLF_ARRAY_FULL_INTEGER
elif dtype ==6:
self.wolftype = WOLF_ARRAY_FULL_SINGLE
elif dtype == 7:
self.wolftype = WOLF_ARRAY_FULL_DOUBLE
else:
logging.error(_('The datatype of the raster is not supported -- {}'.format(dtype)))
logging.error(_('Please convert the raster to a supported datatype - or upgrade the code to support this datatype'))
logging.error(_('See : read_txt_header and import_geotif in wolf_array.py'))
return
# Find the EPSG code
proj = raster.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
if srs.IsProjected:
self.epsg = int(srs.GetAttrValue('AUTHORITY',1))
if "EPSG:" in proj:
try:
idx1 = proj.index("EPSG:")
idx2 = proj.index('"', idx1)
epsg2 = int(proj[idx1+5:idx2])
if epsg2 != self.epsg:
logging.warning(_('Mismatch in EPSG code detection: {} and {}').format(self.epsg, epsg2))
logging.warning(_("Check your file's projection -- You can use gdalinfo command line tool"))
self.epsg = epsg2
logging.warning(_('Using EPSG: {}').format(self.epsg))
except:
logging.warning(_('Could not check EPSG code in the projection string'))
logging.warning(_("Check your file's projection -- You can use gdalinfo command line tool"))
elif filename.endswith('.vrt'):
# Virtual raster
# from osgeo import gdal
raster:gdal.Dataset
raster = gdal.Open(filename)
geotr = raster.GetGeoTransform()
self.dx = geotr[1]
self.dy = abs(geotr[5])
self.origx = geotr[0]
self.origy = geotr[3]
if geotr[5] < 0.:
self.origy = geotr[3] + raster.RasterYSize * geotr[5]
self.nbx = raster.RasterXSize
self.nby = raster.RasterYSize
dtype = raster.GetRasterBand(1).DataType
if dtype == 1:
self.wolftype = WOLF_ARRAY_FULL_INTEGER8
elif dtype in [2,3]:
self.wolftype = WOLF_ARRAY_FULL_INTEGER16
elif dtype in [4,5] :
self.wolftype = WOLF_ARRAY_FULL_INTEGER
elif dtype ==6:
self.wolftype = WOLF_ARRAY_FULL_SINGLE
elif dtype == 7:
self.wolftype = WOLF_ARRAY_FULL_DOUBLE
else:
logging.error(_('The datatype of the raster is not supported -- {}'.format(dtype)))
logging.error(_('Please convert the raster to a supported datatype - or upgrade the code to support this datatype'))
logging.error(_('See : read_txt_header and import_geotif in wolf_array.py'))
return
# Find the EPSG code
proj = raster.GetProjection()
srs = osr.SpatialReference()
srs.ImportFromWkt(proj)
if srs.IsProjected:
self.epsg = int(srs.GetAttrValue('AUTHORITY',1))
if "EPSG:" in proj:
try:
idx1 = proj.index("EPSG:")
idx2 = proj.index('"', idx1)
epsg2 = int(proj[idx1+5:idx2])
if epsg2 != self.epsg:
logging.warning(_('Mismatch in EPSG code detection: {} and {}').format(self.epsg, epsg2))
logging.warning(_("Check your file's projection -- You can use gdalinfo command line tool"))
self.epsg = epsg2
logging.warning(_('Using EPSG: {}').format(self.epsg))
except:
logging.warning(_('Could not check EPSG code in the projection string'))
logging.warning(_("Check your file's projection -- You can use gdalinfo command line tool"))
elif filename.endswith('.npy') and not os.path.exists(filename + '.txt'):
# Il y de fortes chances que cette matrice numpy provienne d'une modélisation GPU
# et donc que les coordonnées et la résolution soient disponibles dans un fichier parameters.json
if (locpath.parent / 'parameters.json').exists():
with open(locpath.parent / 'parameters.json', 'r') as f:
params = json.load(f)
if 'parameters' in params.keys():
if "dx" in params['parameters'].keys() :
self.dx = float(params['parameters']["dx"])
if "dy" in params['parameters'].keys() :
self.dy = float(params['parameters']["dy"])
if "base_coord_x" in params['parameters'].keys() :
self.origx = float(params['parameters']["base_coord_x"])
if "base_coord_y" in params['parameters'].keys() :
self.origy = float(params['parameters']["base_coord_y"])
self.nullvalue = 0 if 'nap.npy' in filename.lower() else 99999.
else:
self.dx = 1.
self.dy = 1.
self.origx = 0.
self.origy = 0.
# Numpy format
with open(filename, 'rb') as f:
version = np.lib.format.read_magic(f)
if version[0] == 1:
shape, fortran, dtype = np.lib.format.read_array_header_1_0(f)
elif version[0] == 2:
shape, fortran, dtype = np.lib.format.read_array_header_2_0(f)
else:
raise ValueError("Unknown numpy version: %s" % version)
self.nbx, self.nby = shape
if dtype == np.float32:
self.wolftype = WOLF_ARRAY_FULL_SINGLE
elif dtype == np.float64:
self.wolftype = WOLF_ARRAY_FULL_DOUBLE
elif dtype == np.int32:
self.wolftype = WOLF_ARRAY_FULL_INTEGER
elif dtype == np.int16:
self.wolftype = WOLF_ARRAY_FULL_INTEGER16
elif dtype == np.uint8:
self.wolftype = WOLF_ARRAY_FULL_UINTEGER8
elif dtype == np.int8:
self.wolftype = WOLF_ARRAY_FULL_INTEGER8
else:
logging.error(_('Unsupported type in numpy file -- Abort loading'))
return
elif filename.endswith('.flt'):
# Fichier .flt
if not os.path.exists(filename[:-4] + '.hdr'):
logging.warning(_('File {} does not exist -- Retry!'.format(filename[:-4] + '.hdr')))
return
f = open(filename[:-4] + '.hdr', 'r')
lines = f.read().splitlines()
f.close()
for curline in lines:
if 'NCOLS' in curline.upper():
tmp = curline.split(' ')
self.nbx = int(tmp[-1])
elif 'NROWS' in curline.upper():
tmp = curline.split(' ')
self.nby = int(tmp[-1])
elif 'XLLCORNER' in curline.upper():
tmp = curline.split(' ')
self.origx = float(tmp[-1])
elif 'YLLCORNER' in curline.upper():
tmp = curline.split(' ')
self.origy = float(tmp[-1])
elif 'ULXMAP' in curline.upper():
tmp = curline.split(' ')
self.origx = float(tmp[-1])
self.flipupd=True
elif 'ULYMAP' in curline.upper():
tmp = curline.split(' ')
self.origy = float(tmp[-1])
self.flipupd=True
elif 'CELLSIZE' in curline.upper():
tmp = curline.split(' ')
self.dx = self.dy = float(tmp[-1])
elif 'XDIM' in curline.upper():
tmp = curline.split(' ')
self.dx = float(tmp[-1])
elif 'YDIM' in curline.upper():
tmp = curline.split(' ')
self.dy = float(tmp[-1])
elif 'NODATA' in curline.upper():
tmp = curline.split(' ')
self.nullvalue = float(tmp[-1])
self.origx -= self.dx*0.5
self.origy += self.dy*0.5
if self.flipupd:
self.origy -= self.dy*float(self.nby)
else:
if not os.path.exists(filename + '.txt'):
logging.info(_('File {} does not exist -- Maybe be a parameter.json exists or retry !'.format(filename + '.txt')))
return
with open(filename + '.txt', 'r') as f:
lines = f.read().splitlines()
tmp = lines[0].split(':')
self.nbx = int(tmp[1])
tmp = lines[1].split(':')
self.nby = int(tmp[1])
tmp = lines[2].split(':')
self.origx = float(tmp[1])
tmp = lines[3].split(':')
self.origy = float(tmp[1])
tmp = lines[4].split(':')
self.dx = float(tmp[1])
tmp = lines[5].split(':')
self.dy = float(tmp[1])
tmp = lines[6].split(':')
self.wolftype = int(tmp[1])
tmp = lines[7].split(':')
self.translx = float(tmp[1])
tmp = lines[8].split(':')
self.transly = float(tmp[1])
decal = 9
if self.wolftype == WOLF_ARRAY_FULL_SINGLE_3D:
tmp = lines[9].split(':')
self.nbz = int(tmp[1])
tmp = lines[10].split(':')
self.origz = float(tmp[1])
tmp = lines[11].split(':')
self.dz = float(tmp[1])
tmp = lines[12].split(':')
self.translz = float(tmp[1])
decal = 13
if self.wolftype in WOLF_ARRAY_MB:
tmp = lines[decal].split(':')
nb_blocks = int(tmp[1])
decal += 1
for i in range(nb_blocks):
curhead = header_wolf()
tmp = lines[decal].split(':')
curhead.nbx = int(tmp[1])
tmp = lines[decal + 1].split(':')
curhead.nby = int(tmp[1])
tmp = lines[decal + 2].split(':')
curhead.origx = float(tmp[1])
tmp = lines[decal + 3].split(':')
curhead.origy = float(tmp[1])
tmp = lines[decal + 4].split(':')
curhead.dx = float(tmp[1])
tmp = lines[decal + 5].split(':')
curhead.dy = float(tmp[1])
decal += 6
curhead.translx = self.translx + self.origx
curhead.transly = self.transly + self.origy
self.head_blocks[getkeyblock(i)] = curhead
@classmethod
[docs]
def read_header(cls, filename:str) -> "header_wolf":
"""
alias for read_txt_header
:param filename: path and filename of the basefile
"""
newhead = cls()
newhead.read_txt_header(filename)
return newhead
[docs]
def write_header(self, filename:str,
wolftype:int,
forceupdate:bool=False):
"""
alias for write_txt_header
"""
self.write_txt_header(filename, wolftype, forceupdate)
[docs]
def write_txt_header(self,
filename:str,
wolftype:int,
forceupdate:bool=False):
"""
Writing the header to a text file
Nullvalue is not written
:param filename: path and filename with '.txt' extension, which will NOT be automatically added
:param wolftype: type of the WOLF_ARRAY_* array
:param forceupdate: if True, the file is rewritten even if it already exists
"""
assert wolftype in [WOLF_ARRAY_CSR_DOUBLE, WOLF_ARRAY_FULL_SINGLE, WOLF_ARRAY_FULL_DOUBLE, WOLF_ARRAY_SYM_DOUBLE, WOLF_ARRAY_FULL_LOGICAL,
WOLF_ARRAY_CSR_DOUBLE, WOLF_ARRAY_FULL_INTEGER, WOLF_ARRAY_FULL_SINGLE_3D, WOLF_ARRAY_FULL_INTEGER8, WOLF_ARRAY_FULL_UINTEGER8,
WOLF_ARRAY_MB_SINGLE, WOLF_ARRAY_MB_INTEGER, WOLF_ARRAY_FULL_INTEGER16, WOLF_ARRAY_MNAP_INTEGER, WOLF_ARRAY_FULL_INTEGER16_2], _('The type of array is not correct')
if not os.path.exists(filename) or forceupdate:
with open(filename,'w') as f:
""" Ecriture de l'en-tête de Wolf array """
f.write('NbX :\t{0}\n'.format(str(self.nbx)))
f.write('NbY :\t{0}\n'.format(str(self.nby)))
f.write('OrigX :\t{0}\n'.format(str(self.origx)))
f.write('OrigY :\t{0}\n'.format(str(self.origy)))
f.write('DX :\t{0}\n'.format(str(self.dx)))
f.write('DY :\t{0}\n'.format(str(self.dy)))
f.write('TypeEnregistrement :\t{0}\n'.format(str(wolftype)))
f.write('TranslX :\t{0}\n'.format(str(self.translx)))
f.write('TranslY :\t{0}\n'.format(str(self.transly)))
if wolftype == WOLF_ARRAY_FULL_SINGLE_3D:
f.write('NbZ :\t{0}\n'.format(str(self.nbz)))
f.write('OrigZ :\t{0}\n'.format(str(self.origz)))
f.write('DZ :\t{0}\n'.format(str(self.dz)))
f.write('TranslZ :\t{0}\n'.format(str(self.translz)))
if wolftype in WOLF_ARRAY_MB:
f.write('Nb Blocs :\t{0}\n'.format(str(self.nb_blocks)))
for i in range(self.nb_blocks):
curhead = self.head_blocks[getkeyblock(i)]
f.write('NbX :\t{0}\n'.format(str(curhead.nbx)))
f.write('NbY :\t{0}\n'.format(str(curhead.nby)))
f.write('OrigX :\t{0}\n'.format(str(curhead.origx)))
f.write('OrigY :\t{0}\n'.format(str(curhead.origy)))
f.write('DX :\t{0}\n'.format(str(curhead.dx)))
f.write('DY :\t{0}\n'.format(str(curhead.dy)))
[docs]
def is_like(self, other:"header_wolf", check_mb:bool=False) -> bool:
"""
Comparison of two headers
:param other: other header to compare
:param check_mb: if True, the comparison is done on the blocks too
The nullvalue is not taken into account
"""
test = True
test &= self.origx == other.origx
test &= self.origy == other.origy
test &= self.origz == other.origz
test &= self.translx == other.translx
test &= self.transly == other.transly
test &= self.translz == other.translz
test &= self.dx == other.dx
test &= self.dy == other.dy
test &= self.dz == other.dz
test &= self.nbx == other.nbx
test &= self.nby == other.nby
test &= self.nbz == other.nbz
test &= self.nbdims == other.nbdims
if check_mb:
test &= self.nb_blocks == other.nb_blocks
for block1, block2 in zip(self.head_blocks.values(), other.head_blocks.values()):
test &= block1.is_like(block2)
return test
[docs]
def align2grid(self, x1:float, y1:float, eps:float=0.0001) -> tuple[float,float]:
"""
Align coordinates to nearest grid point
where the grid is defined by the borders of the array.
"""
if x1-self.origx < 0:
x2 = np.round((x1 - self.origx + eps) / self.dx) * self.dx + self.origx
else:
x2 = np.round((x1 - self.origx - eps) / self.dx) * self.dx + self.origx
if y1-self.origy < 0:
y2 = np.round((y1 - self.origy + eps) / self.dy) * self.dy + self.origy
else:
y2 = np.round((y1 - self.origy - eps) / self.dy) * self.dy + self.origy
return x2, y2
[docs]
def _rasterize_segment(self,
x1:float, y1:float,
x2:float, y2:float,
xstart:float=None,
ystart:float=None,
n_packet:int = 10) -> list[list[float]]:
"""
Rasterize a segment according to the grid
where the grid is defined by the borders of the array.
:param x1: x coordinate of the first point
:param y1: y coordinate of the first point
:param x2: x coordinate of the second point
:param y2: y coordinate of the second point
:param xstart: x coordinate of the starting point
:param ystart: y coordinate of the starting point
:return: numpy array of the rasterized segment
"""
if xstart is None and ystart is None:
xstart, ystart = self.align2grid(x1, y1)
x2, y2 = self.align2grid(x2, y2)
points=[]
points.append([xstart, ystart])
length = 99999.
prec = min(self.dx, self.dy)
direction = np.array([x2-xstart, y2-ystart])
length = np.linalg.norm(direction)
if length == 0.:
logging.debug(_('Segment of zero length -- return starting point only'))
return points
direction /= length
while length >= prec:
if np.abs(direction[0])>= np.abs(direction[1]):
if length > n_packet*prec:
# I we are far from the end of the segment, we can use the n_packet to rasterize
# this will be used to avoid to many straight lines
nx = np.abs(int(direction[0] * (n_packet-1))) # number of points in x direction
ny = np.abs(int(direction[1] * (n_packet-1))) # number of points in y direction
n_common = min(nx, ny) # number of common points
for i in range(n_common):
xstart += self.dx * np.sign(direction[0])
points.append([xstart, ystart])
ystart += self.dy * np.sign(direction[1])
points.append([xstart, ystart])
for i in range(nx - n_common):
xstart += self.dx * np.sign(direction[0])
points.append([xstart, ystart])
for j in range(ny - n_common):
ystart += self.dy * np.sign(direction[1])
points.append([xstart, ystart])
else:
xstart += self.dx * np.sign(direction[0])
points.append([xstart, ystart])
else:
if length > n_packet*prec:
nx = np.abs(int(direction[0] * (n_packet-1)))
ny = np.abs(int(direction[1] * (n_packet-1)))
n_common = min(nx, ny)
for j in range(n_common):
ystart += self.dy * np.sign(direction[1])
points.append([xstart, ystart])
xstart += self.dx * np.sign(direction[0])
points.append([xstart, ystart])
for j in range(ny - n_common):
ystart += self.dy * np.sign(direction[1])
points.append([xstart, ystart])
for i in range(nx - n_common):
xstart += self.dx * np.sign(direction[0])
points.append([xstart, ystart])
else:
ystart += self.dy * np.sign(direction[1])
points.append([xstart, ystart])
direction = np.array([x2-xstart, y2-ystart])
length = np.linalg.norm(direction)
if length > 0.:
direction /= length
return points
[docs]
def rasterize_vector_along_grid(self, vector2raster:"vector", outformat=None) -> Union[np.ndarray,"vector"]:
"""
Rasterize a vector according to the grid
:param vector2raster: vector to rasterize
:param outformat: output format (np.ndarray or vector)
"""
from ..PyVertexvectors import vector as _vector
if outformat is None:
outformat = _vector
assert outformat in [np.ndarray, _vector], _('outformat must be np.ndarray or vector')
# get the vertices of the vector
xy = vector2raster.asnparray().tolist()
# rasterize the vector
rasterized = []
rasterized += self._rasterize_segment(xy[0][0], xy[0][1], xy[1][0], xy[1][1])
for i in range(1, len(xy)-1):
out = self._rasterize_segment(xy[i][0], xy[i][1],
xy[i+1][0], xy[i+1][1],
rasterized[-1][0], rasterized[-1][1])
if len(out) > 1:
rasterized += out[1:]
# get the indices of the rasterized vector
xy = np.array(rasterized)
if outformat is np.ndarray:
return xy
elif outformat is _vector:
#create new vector
newvector = _vector()
newvector.add_vertices_from_array(xy)
return newvector
[docs]
def rasterize_zone_along_grid(self, zone2raster:"zone", outformat=None) -> Union[np.ndarray,"zone"]:
""" Rasterize a zone according to the grid
:param zone2raster: zone to rasterize
:param outformat: output format (np.ndarray or zone)
"""
from ..PyVertexvectors import zone as _zone
if outformat is None:
outformat = _zone
assert outformat in [np.ndarray, _zone], _('outformat must be np.ndarray or zone')
if outformat is _zone:
out_vec = None # will default to vector in rasterize_vector_along_grid
else:
out_vec = np.ndarray
new_vecs = []
for vec in zone2raster.myvectors:
new_vecs.append(self.rasterize_vector_along_grid(vec, outformat=out_vec))
if outformat is np.ndarray:
return np.vstack(new_vecs)
elif outformat is _zone:
# create new zone
newzone = _zone()
for vec, oldvec in zip(new_vecs, zone2raster.myvectors):
vec.myname = oldvec.myname # keep the original name
newzone.add_vector(vec, forceparent=True)
return newzone
[docs]
def rasterize_vector(self, vector2raster:"vector", outformat=None) -> Union[np.ndarray,"vector"]:
""" DEPRECATED since 2.2.8 -- use rasterize_vector_along_grid instead.
Will be removed in 2.3.0
"""
logging.warning(_('rasterize_vector is deprecated since 2.2.8 -- use rasterize_vector_along_grid instead'))
return self.rasterize_vector_along_grid(vector2raster, outformat=outformat)
[docs]
def get_xy_infootprint_vect(self, myvect: "vector | Polygon", eps:float = 0.) -> tuple[np.ndarray,np.ndarray]:
"""
Return the coordinates and the indices of the cells in the footprint of a vector.
Coordinates are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint.
Indices are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint. -> See get_ij_infootprint
:param myvect: target vector or Shapely Polygon
:param eps: epsilon to avoid rounding errors
:return: tuple of two numpy arrays - (coordinates, indices)
"""
myptsij = self.get_ij_infootprint_vect(myvect, eps=eps)
mypts=np.asarray(myptsij.copy(),dtype=np.float64)
mypts[:,0] = (mypts[:,0]+.5)*self.dx +self.origx +self.translx
mypts[:,1] = (mypts[:,1]+.5)*self.dy +self.origy +self.transly
return mypts, myptsij
[docs]
def get_xy_infootprint(self, myvect: "vector | Polygon", eps:float = 0.) -> tuple[np.ndarray,np.ndarray]:
"""
Return the coordinates and the indices of the cells in the footprint of a vector.
Coordinates are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint.
Indices are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint. -> See get_ij_infootprint
Main principle:
- get the indices with 'get_ij_infootprint'
- then convert them to coordinates.
:param myvect: target vector or Shapely Polygon
:param eps: epsilon to avoid rounding errors
:return: tuple of two numpy arrays - (coordinates, indices)
"""
return self.get_xy_infootprint_vect(myvect, eps=eps)
[docs]
def get_ij_infootprint_vect(self, myvect: "vector | Polygon", eps:float = 0.) -> np.ndarray:
"""
Return the indices of the cells in the footprint of a vector.
Main principle:
- get the bounding box of the vector
- convert the bounding box to indices
- limit indices to the array size
- create a meshgrid of indices
Indices are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint.
:param myvect: target vector or Shapely Polygon
:param eps: epsilon to avoid rounding errors
:return: numpy array of indices
"""
from ..PyVertexvectors import vector as _vector
if isinstance(myvect, Polygon):
xmin, ymin, xmax, ymax = myvect.bounds
elif isinstance(myvect, _vector):
xmin, ymin, xmax, ymax = myvect.xmin, myvect.ymin, myvect.xmax, myvect.ymax
else:
logging.error(_('The object must be a vector or a Polygon'))
return np.array([])
i1, j1 = self.get_ij_from_xy(xmin+eps, ymin+eps)
i2, j2 = self.get_ij_from_xy(xmax-eps, ymax-eps)
i1 = max(i1,0) # FIXME Why ??? How could i,j be negative ? --> because this function can be called with a vector that is not in the array (e.g. a vector defined by clicks in the UI)
j1 = max(j1,0)
i2 = min(i2,self.nbx-1)
j2 = min(j2,self.nby-1)
xv,yv = np.meshgrid(np.arange(i1,i2+1),np.arange(j1,j2+1))
mypts = np.hstack((xv.flatten()[:,np.newaxis],yv.flatten()[:,np.newaxis]))
return mypts
[docs]
def get_ij_infootprint(self, myvect: "vector | Polygon", eps:float = 0.) -> np.ndarray:
"""
Return the indices of the cells in the footprint of a vector
Main principle:
- get the bounding box of the vector
- convert the bounding box to indices
- limit indices to the array size
- create a meshgrid of indices
Indices are stored in a numpy array of shape (n,2) where n is the number of cells in the footprint.
:param myvect: target vector or Shapely Polygon
:param eps: epsilon to avoid rounding errors
:return: numpy array of indices
"""
return self.get_ij_infootprint_vect(myvect, eps=eps)
[docs]
def convert_xy2ij_np(self,xy):
"""
Convert XY coordinates to IJ indices **(0-based)** with Numpy without any check/options
:param xy: = numpy array of shape (n,2) with XY coordinates
"""
return np.asarray((xy[:,0]-self.origx -self.translx)/self.dx-.5,dtype=np.int32), \
np.asarray((xy[:,1]-self.origy -self.transly)/self.dy-.5,dtype=np.int32)
[docs]
def convert_ij2xy_np(self,ij):
"""
Convert IJ indices **(0-based)** to XY coordinates with Numpy without any check/options
:param ij: = numpy array of shape (n,2) with IJ indices
"""
return np.asarray((ij[:,0]+.5)*self.dx+self.origx +self.translx ,dtype=np.float64), \
np.asarray((ij[:,1]+.5)*self.dy+self.origy +self.transly ,dtype=np.float64)
[docs]
def is_in_footprint(self, x:float, y:float, eps:float=0.) -> bool:
"""
Check if a point is in the footprint of the array
:param x: x coordinate
:param y: y coordinate
:param eps: epsilon to avoid rounding errors
:return: True if the point is in the footprint, False otherwise
"""
return (x >= self.origx + self.translx - eps) and (x <= self.origx + self.translx + self.nbx * self.dx + eps) and \
(y >= self.origy + self.transly - eps) and (y <= self.origy + self.transly + self.nby * self.dy + eps)
@classmethod
[docs]
def make(self, orig: tuple[float, float] = (0.0, 0.0), nb: tuple[int, int]=(0,0), d:tuple[float, float]=(1.0,1.0)):
wh = header_wolf()
wh.nbx, wh.nby = nb
wh.origx, wh.origy = orig
wh.dx, wh.dy = d
return wh
@classmethod
[docs]
def make_from_xybounds_grid(self,
xbounds:tuple[float,float], ybounds:tuple[float,float],
d:tuple[float,float]=(1.0,1.0), grid_to_align:'header_wolf' = None):
""" Create a header_wolf from the bounds in x and y and the resolution in x and y.
If grid_size is provided, the origin will be adjusted to fit the grid size.
:param xbounds: tuple of two floats - (xmin, xmax)
:param ybounds: tuple of two floats - (ymin, ymax)
:param d: tuple of two floats - (dx, dy)
:param grid_to_align: header_wolf - grid to adjust the origin
:return: header_wolf instance
"""
wh = header_wolf()
wh.dx, wh.dy = d
# The origin are at the lower left corner.
# So we need to adjust the bounds accordingly the grid size.
if grid_to_align is not None:
wh.origx = (xbounds[0] - ((xbounds[0] - grid_to_align.origx) % grid_to_align.dx))
wh.origy = (ybounds[0] - ((ybounds[0] - grid_to_align.origy) % grid_to_align.dy))
else:
wh.origx = xbounds[0]
wh.origy = ybounds[0]
wh.nbx = int(np.ceil( (xbounds[1]-wh.origx) / wh.dx ))
wh.nby = int(np.ceil( (ybounds[1]-wh.origy) / wh.dy ))
return wh
@classmethod
[docs]
def make_from_header(self, other:"header_wolf"):
wh = header_wolf()
wh.nbx, wh.nby = other.nbx, other.nby
wh.origx, wh.origy = other.origx, other.origy
wh.dx, wh.dy = other.dx, other.dy
return wh
[docs]
def copy(self) -> "header_wolf":
""" Return a copy of the header_wolf instance """
return header_wolf.make_from_header(self)
@classmethod
[docs]
def from_slices(self, x_slice:slice, y_slice:slice, d:tuple[float, float]=(1.0,1.0)):
wh = header_wolf()
wh.dx, wh.dy = d
wh.nbx, wh.nby = x_slice.stop - x_slice.start, y_slice.stop - y_slice.start
wh.origx, wh.origy = x_slice.start * wh.dx, y_slice.start * wh.dy
return wh
[docs]
def to_slices(self) -> tuple[slice, slice]:
x_slice = slice(int(self.origx / self.dx), int(self.origx / self.dx) + self.nbx)
y_slice = slice(int(self.origy / self.dy), int(self.origy / self.dy) + self.nby)
return x_slice, y_slice
[docs]
def to_dict(self, type_to_save: Literal["header_wolf", "bounds"] = "header_wolf") -> dict:
if type_to_save == "header_wolf":
d= {
"header_type": "header_wolf",
"dx": self.dx,
"dy": self.dy,
"origx": self.origx,
"origy": self.origy,
"nbx": self.nbx,
"nby": self.nby
}
elif type_to_save == "bounds":
d= {
"header_type": "bounds",
"minx": self.origx,
"miny": self.origy,
"maxx": self.origx + self.dx * self.nbx,
"maxy": self.origy + self.dy * self.nby,
}
else:
raise ValueError(f"Unknown type_to_save '{type_to_save}'. Supported values are 'header_wolf' and 'bounds'.")
return d
@classmethod
[docs]
def from_dict(cls, d:dict, dxdy:tuple[float,float]=None, grid:"header_wolf"=None) -> "header_wolf":
hdr_type = d.get("header_type", "header_wolf")
newhdr = header_wolf()
if hdr_type == "header_wolf":
return header_wolf.make(
orig=(d["origx"], d["origy"]),
nb=(d["nbx"], d["nby"]),
d=(d["dx"], d["dy"])
)
elif hdr_type == "bounds":
assert dxdy is not None or grid is not None, _('dxdy must be provided when loading from bounds')
if dxdy is not None:
if grid is None:
grid = header_wolf.make(orig = (d["minx"], d["miny"]),
nb = (1,1),
d = dxdy)
return header_wolf.make_from_xybounds_grid(xbounds=(d["minx"], d["maxx"]),
ybounds=(d["miny"], d["maxy"]),
d=dxdy,
grid_to_align=grid)
if grid is not None:
return header_wolf.make_from_xybounds_grid(
xbounds=(d["minx"], d["maxx"]),
ybounds=(d["miny"], d["maxy"]),
d=dxdy,
grid_to_align=grid
)
elif grid is not None:
return header_wolf.make_from_xybounds_grid(
xbounds=(d["minx"], d["maxx"]),
ybounds=(d["miny"], d["maxy"]),
d=(grid.dx, grid.dy),
grid_to_align=grid
)