{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Entities\n",
    "\n",
    "This section introduces the different entities that can be created and stored in the `geoh5` file format.\n",
    "\n",
    "![entities](./images/entities.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Groups\n",
    "\n",
    "`Groups` are effectively containers for other entities, such as ``Objects`` (Points, Curve, Surface, etc.) and other `Groups`. Groups are used to establish `parent-child` relationships and to store information about a collection of entities.  \n",
    "\n",
    "### RootGroup\n",
    "By default, the parent of any new `Entity` is the workspace ``RootGroup``.  It is the only entity in the ``Workspace`` without a parent. Users rarely have to interect with the ``Root`` group as it is mainly used to maintain the overall project hierarchy.\n",
    "\n",
    "![Root](./images/root.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ContainerGroup\n",
    "\n",
    "\n",
    "A ``ContainerGroup`` can easily be added to the workspace and can be assigned a `name` and `description`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.groups import ContainerGroup\n",
    "from geoh5py.workspace import Workspace\n",
    "import numpy as np\n",
    "\n",
    "# Create a blank project\n",
    "workspace = Workspace(\"my_project.geoh5\")\n",
    "\n",
    "# Add a group\n",
    "group = ContainerGroup.create(workspace, name='myGroup')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At creation, `\"myGroup\"` is written to the project ``geoh5`` file and visible in the Analyst project tree.\n",
    "\n",
    "![Groups](./images/groups.png)\n",
    "\n",
    "Any entity can be accessed by its `name` or `uid` (unique identifier):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(group.uid)\n",
    "print(workspace.get_entity(\"myGroup\")[0] == workspace.get_entity(group.uid)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objects\n",
    "\n",
    "The `geoh5` format enables storing a wide variety of ``Object`` entities that can be displayed in 3D. \n",
    "This section describes the collection of ``Objects`` entities currently supported by `geoh5py`.\n",
    "\n",
    "![Gobjects](./images/objects.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Points\n",
    "\n",
    "The ``Points`` object consists of a list of `vertices` that define the location of actual data in 3D space. As for all other `Objects`, it can be created from an array of 3D coordinates and added to any group as follow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import Points\n",
    "\n",
    "# Generate a numpy array of xyz locations\n",
    "n = 100\n",
    "radius, theta = np.arange(n), np.linspace(0, np.pi*8, n)\n",
    "\n",
    "x, y = radius * np.cos(theta), radius * np.sin(theta)\n",
    "z = (x**2. + y**2.)**0.5\n",
    "xyz = np.c_[x.ravel(), y.ravel(), z.ravel()] # Form a 2D array\n",
    "\n",
    "# Create the Point object\n",
    "points = Points.create(\n",
    "    workspace,       # The target Workspace\n",
    "    vertices=xyz     # Set vertices\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![points](./images/points.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Curve\n",
    "\n",
    "The ``Curve`` object, also known as a polyline, is often used to define contours, survey lines or geological contacts. It is a sub-class of the ``Points`` object with the added `cells` property, that defines the line segments connecting its `vertices`. By default, all vertices are connected sequentially following the order of the input `vertices`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import Curve\n",
    "\n",
    "# Create the Curve object\n",
    "curve = Curve.create(\n",
    "    workspace,       # The target Workspace\n",
    "    vertices=xyz\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, the `cells` property can be modified, either directly or by assigning `parts` identification to each `vertices`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the curve into two parts\n",
    "part_id = np.ones(n, dtype=\"int32\")\n",
    "part_id[:75] = 2\n",
    "\n",
    "# Assign the part\n",
    "curve.parts = part_id"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![line](./images/line.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drillhole\n",
    "\n",
    "\n",
    "``Drillhole`` objects are different from other objects as their 3D geometry is defined by the `collar` and `surveys` attributes. As for version `geoh5 v2.0`, the drillholes require a `DrillholeGroup` entity to store the geometry and data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.groups import DrillholeGroup\n",
    "from geoh5py.objects import Drillhole\n",
    "\n",
    "dh_group = DrillholeGroup.create(workspace)\n",
    "\n",
    "# Create a simple well\n",
    "total_depth = 100\n",
    "dist = np.linspace(0, total_depth, 10)\n",
    "azm = np.ones_like(dist) * 45.\n",
    "dip = np.linspace(-89, -75, dist.shape[0])\n",
    "collar = np.r_[0., 10., 10]\n",
    "\n",
    "well = Drillhole.create(\n",
    "    workspace, collar=collar, surveys=np.c_[dist, azm, dip], name=\"Drillhole\", parent=dh_group \n",
    ")\n",
    "\n",
    "print(well.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![drillhole](./images/drillhole.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Surface\n",
    "\n",
    "\n",
    "The ``Surface`` object is also described by `vertices` and `cells` that form a net of triangles. If omitted on creation, the `cells` property is calculated using a 2D `scipy.spatial.Delaunay` triangulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import Surface\n",
    "from scipy.spatial import Delaunay\n",
    "\n",
    "# Create a triangulated surface from points\n",
    "surf_2D = Delaunay(xyz[:, :2])\n",
    "\n",
    "# Create the Surface object\n",
    "surface = Surface.create(\n",
    "    workspace, \n",
    "    vertices=points.vertices, # Add vertices\n",
    "    cells=surf_2D.simplices\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![surface](./images/surface.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GeoImage\n",
    "\n",
    "The `GeoImage` object handles raster data, either single or 3-band images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import GeoImage\n",
    "\n",
    "geoimage = GeoImage.create(workspace)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Image values can be assigned to the object from either a 2D `numpy.ndarray` for single band (gray):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geoimage.image = np.random.randn(128, 128)\n",
    "display(geoimage.image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or as 3D `numpy.ndarray` for 3-band RGB image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geoimage.image = np.random.randn(128, 128, 3)\n",
    "display(geoimage.image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or directly from file (png, jpeg, tiff)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "geoimage.image = \"./images/flin_flin_geology.jpg\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `PIL.Image` object gets exposed to the user, which can be used for common raster manipulation (rotation, filtering, etc). The modified raster is stored back on file as a blob (`bytes`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(geoimage.image)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Geo-referencing\n",
    "By default, the `GeoImage` object will be displayed at the origin (xy-plane) with dimensions equal to the pixel count. The utility function [GeoImage.georeference](../api/geoh5py.objects.rst#geoh5py.objects.geo_image.GeoImage.georeference) lets users geo-reference the image in 3D space based on at least three (3) input reference points (pixels) with associated world coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixels = [\n",
    "    [18, 73],\n",
    "    [757, 1014],\n",
    "    [18, 1014],\n",
    "\n",
    "]\n",
    "coords = [\n",
    "    [311005, 6065252, 0],\n",
    "    [320001, 6076748, 0],\n",
    "    [311005, 6076748, 0]\n",
    "]\n",
    "\n",
    "geoimage.georeference(pixels, coords) \n",
    "print(geoimage.vertices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Grid2D\n",
    "\n",
    "\n",
    "The `Grid2D` object defines a regular grid of `cells` often used to display model sections or to compute data derivatives.\n",
    "A `Grid2D` can be oriented in 3D space using the `origin`, `rotation` and `dip` parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import Grid2D\n",
    "\n",
    "# Create the Surface object\n",
    "grid = Grid2D.create(\n",
    "    workspace,\n",
    "    origin = [25, -75, 50],\n",
    "    u_cell_size = 2.5,\n",
    "    v_cell_size = 2.5,\n",
    "    u_count = 64,\n",
    "    v_count = 16,\n",
    "    rotation = 90.0,\n",
    "    dip = 45.0,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![grid2d](./images/grid2d.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DrapeModel\n",
    "\n",
    "The `DrapeModel` object defines an array of vertical (curtain) cells draped below a curved trace. \n",
    "\n",
    "- The `prisms` attribute defines the elevation and position of the uppermost cell faces.\n",
    "- The `layers` attribute defines the bottom face elevation of cells. \n",
    "\n",
    "In the example below we create a simple `DrapeModel` object along a sinusoidal path with equal number of layers below every station. Variable number of layers per prism is also supported.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import DrapeModel\n",
    "\n",
    "# Define a simple trace\n",
    "n_columns = 32\n",
    "n_layers = 8\n",
    "\n",
    "x = np.linspace(0, np.pi, n_columns)\n",
    "y = np.cos(x)\n",
    "z = np.linspace(0, 1, n_columns)\n",
    "\n",
    "# Count the index and number of values per columns\n",
    "layer_count = np.ones(n_columns) * n_layers\n",
    "prisms = np.c_[\n",
    "    x, \n",
    "    y, \n",
    "    z, \n",
    "    np.cumsum(layer_count) - n_layers,\n",
    "    layer_count\n",
    "]\n",
    "\n",
    "# Define the index and elevation of draped cells\n",
    "k_index, i_index = np.meshgrid(np.arange(n_layers), np.arange(n_columns))\n",
    "z_elevation = z[i_index] - np.linspace(0.5, 2, n_layers)\n",
    "\n",
    "layers = np.c_[\n",
    "    i_index.flatten(), k_index.flatten(), z_elevation.flatten()\n",
    "]\n",
    "\n",
    "# Create the object\n",
    "drape = DrapeModel.create(workspace, prisms=prisms, layers=layers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![drapemodel](./images/drapemodel.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BlockModel\n",
    "\n",
    "\n",
    "The ``BlockModel`` object defines a rectilinear grid of cells, also known as a tensor mesh. The `cells` center position is determined by `cell_delimiters` (offsets) along perpendicular axes (`u`, `v`, `z`) and relative to the `origin`. ``BlockModel`` can be oriented horizontally by controlling the `rotation` parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import BlockModel\n",
    "\n",
    "# Create the Surface object\n",
    "blockmodel = BlockModel.create(\n",
    "    workspace,\n",
    "    origin = [25, -100, 50],\n",
    "    u_cell_delimiters=np.cumsum(np.ones(16) * 5), # Offsets along u\n",
    "    v_cell_delimiters=np.cumsum(np.ones(32) * 5), # Offsets along v\n",
    "    z_cell_delimiters=np.cumsum(np.ones(16) * -2.5),  # Offsets along z (down)\n",
    "    rotation = 30.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![blockmodel](./images/blockmodel.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Octree\n",
    "\n",
    "\n",
    "The ``Octree`` object is type of 3D grid that uses a tree structure to define `cells`. Each cell can be subdivided it into eight octants allowing for a more efficient local refinement of the mesh. The ``Octree`` object can also be oriented horizontally by controlling the `rotation` parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoh5py.objects import Octree\n",
    "\n",
    "octree = Octree.create(\n",
    "        workspace,\n",
    "        origin=[25, -100, 50],\n",
    "        u_count=16,      # Number of cells in power 2\n",
    "        v_count=32,      \n",
    "        w_count=16,\n",
    "        u_cell_size=5.0, # Base cell size (highest octree level)\n",
    "        v_cell_size=5.0,\n",
    "        w_cell_size=2.5, # Offsets along z (down)\n",
    "        rotation=30,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the octree mesh will be refined at the lowest level possible along each axes.\n",
    "\n",
    "![octree](./images/octree.png){width=\"50%\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "workspace.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}