diff --git a/README.md b/README.md
index 04798cc..ec1d67e 100644
--- a/README.md
+++ b/README.md
@@ -61,22 +61,26 @@ Usage: ctb-tile [options] GDAL_DATASOURCE
 
 Options:
 
-  -V, --version                 output program version
-  -h, --help                    output help information
-  -o, --output-dir <dir>        specify the output directory for the tiles (defaults to working directory)
-  -f, --output-format <format>  specify the output format for the tiles. This is either `Terrain` (the default) or any format listed by `gdalinfo --formats`
-  -p, --profile <profile>       specify the TMS profile for the tiles. This is either `geodetic` (the default) or `mercator`
-  -c, --thread-count <count>    specify the number of threads to use for tile generation. On multicore machines this defaults to the number of CPUs
-  -t, --tile-size <size>        specify the size of the tiles in pixels. This defaults to 65 for terrain tiles and 256 for other GDAL formats
-  -s, --start-zoom <zoom>       specify the zoom level to start at. This should be greater than the end zoom level
-  -e, --end-zoom <zoom>         specify the zoom level to end at. This should be less than the start zoom level and >= 0
-  -r, --resampling-method <algorithm> specify the raster resampling algorithm.  One of: nearest; bilinear; cubic; cubicspline; lanczos; average; mode; max; min; med; q1; q3. Defaults to average.
-  -n, --creation-option <option> specify a GDAL creation option for the output dataset in the form NAME=VALUE. Can be specified multiple times. Not valid for Terrain tiles.
-  -z, --error-threshold <threshold> specify the error threshold in pixel units for transformation approximation. Larger values should mean faster transforms. Defaults to 0.125
-  -m, --warp-memory <bytes>     The memory limit in bytes used for warp operations. Higher settings should be faster. Defaults to a conservative GDAL internal setting.
-  -R, --resume                  Do not overwrite existing files
-  -q, --quiet                   only output errors
-  -v, --verbose                 be more noisy
+  -V, --version                       output program version
+  -h, --help                          output help information
+  -o --output-dir <dir>               specify the output directory for the tiles (defaults to working directory)
+  -f --output-format <format>         specify the output format for the tiles. This is either `Terrain` (the default), `Mesh` (Chunked LOD mesh), or any format listed by `gdalinfo --formats`
+  -p --profile <profile>              specify the TMS profile for the tiles. This is either `geodetic` (the default) or `mercator`
+  -c --thread-count <count>           specify the number of threads to use for tile generation. On multicore machines this defaults to the number of CPUs
+  -t --tile-size <size>               specify the size of the tiles in pixels. This defaults to 65 for terrain tiles and 256 for other GDAL formats
+  -s --start-zoom <zoom>              specify the zoom level to start at. This should be greater than the end zoom level
+  -e --end-zoom <zoom>                specify the zoom level to end at. This should be less than the start zoom level and >= 0
+  -r --resampling-method <algorithm>  specify the raster resampling algorithm.  One of: nearest; bilinear; cubic; cubicspline; lanczos; average; mode; max; min; med; q1; q3. Defaults to average.
+  -n --creation-option <option>       specify a GDAL creation option for the output dataset in the form NAME=VALUE. Can be specified multiple times. Not valid for Terrain tiles.
+  -z --error-threshold <threshold>    specify the error threshold in pixel units for transformation approximation. Larger values should mean faster transforms. Defaults to 0.125
+  -m --warp-memory <bytes>            specify the memory limit in bytes used for warp operations. Higher settings should be faster. Defaults to a conservative GDAL internal setting.
+  -R --resume                         flag do not overwrite existing files
+  -g --mesh-qfactor <factor>          specify the factor to multiply the estimated geometric error to convert heightmaps to irregular meshes. Larger values should mean minor quality. Defaults to 1.0
+  -l --layer                          flag only outputs the layer.json metadata file
+  -C --cesium-friendly                flag forces the creation of missing root tiles to be CesiumJS-friendly
+  -N --vertex-normals                 flag writes 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting, only for `Mesh` format
+  -q --quiet                          flag outputs only errors
+  -v --verbose                        flag outputs more noisy
 ```
 
 #### Recommendations
diff --git a/src/BoundingSphere.hpp b/src/BoundingSphere.hpp
new file mode 100644
index 0000000..a082f91
--- /dev/null
+++ b/src/BoundingSphere.hpp
@@ -0,0 +1,184 @@
+#ifndef BBSPHERE_HPP
+#define BBSPHERE_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file BoundingSphere.hpp
+ * @brief This declares and defines the `BoundingSphere` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include <vector>
+#include <limits>
+
+#include "Coordinate3D.hpp"
+#include "types.hpp"
+
+namespace ctb {
+  template <class T> class BoundingSphere;
+  template <class T> class BoundingBox;
+}
+
+/// A spherical bounding region which is defined by a center point and a radius
+template <class T>
+class ctb::BoundingSphere {
+public:
+  Coordinate3D<T> center; ///< The center of the BoundingSphere
+  double radius; ///< The radius of the BoundingSphere
+
+  /// Create an empty BoundingSphere
+  BoundingSphere() {
+  }
+  /// Create a BoundingSphere from the specified point stream
+  BoundingSphere(const std::vector<Coordinate3D<T>> &points) {
+    fromPoints(points);
+  }
+
+  /// Calculate the center and radius from the specified point stream
+  /// Based on Ritter's algorithm
+  void fromPoints(const std::vector<Coordinate3D<T>> &points) {
+    const T MAX =  std::numeric_limits<T>::infinity();
+    const T MIN = -std::numeric_limits<T>::infinity();
+
+    Coordinate3D<T> minPointX(MAX, MAX, MAX);
+    Coordinate3D<T> minPointY(MAX, MAX, MAX);
+    Coordinate3D<T> minPointZ(MAX, MAX, MAX);
+    Coordinate3D<T> maxPointX(MIN, MIN, MIN);
+    Coordinate3D<T> maxPointY(MIN, MIN, MIN);
+    Coordinate3D<T> maxPointZ(MIN, MIN, MIN);
+
+    // Store the points containing the smallest and largest component
+    // Used for the naive approach
+    for (int i = 0, icount = points.size(); i < icount; i++) {
+      const Coordinate3D<T> &point = points[i];
+
+      if (point.x < minPointX.x) minPointX = point;
+      if (point.y < minPointY.y) minPointY = point;
+      if (point.z < minPointZ.z) minPointZ = point;
+      if (point.x > maxPointX.x) maxPointX = point;
+      if (point.y > maxPointY.y) maxPointY = point;
+      if (point.z > maxPointZ.z) maxPointZ = point;
+    }
+
+    // Squared distance between each component min and max
+    T xSpan = (maxPointX - minPointX).magnitudeSquared();
+    T ySpan = (maxPointY - minPointY).magnitudeSquared();
+    T zSpan = (maxPointZ - minPointZ).magnitudeSquared();
+
+    Coordinate3D<T> diameter1 = minPointX;
+    Coordinate3D<T> diameter2 = maxPointX;
+    T maxSpan = xSpan;
+    if (ySpan > maxSpan) {
+      diameter1 = minPointY;
+      diameter2 = maxPointY;
+      maxSpan = ySpan;
+    }
+    if (zSpan > maxSpan) {
+      diameter1 = minPointZ;
+      diameter2 = maxPointZ;
+      maxSpan = zSpan;
+    }
+
+    Coordinate3D<T> ritterCenter = Coordinate3D<T>(
+      (diameter1.x + diameter2.x) * 0.5,
+      (diameter1.y + diameter2.y) * 0.5,
+      (diameter1.z + diameter2.z) * 0.5
+    );
+    T radiusSquared = (diameter2 - ritterCenter).magnitudeSquared();
+    T ritterRadius = std::sqrt(radiusSquared);
+
+    // Initial center and radius (naive) get min and max box
+    Coordinate3D<T> minBoxPt(minPointX.x, minPointY.y, minPointZ.z);
+    Coordinate3D<T> maxBoxPt(maxPointX.x, maxPointY.y, maxPointZ.z);
+    Coordinate3D<T> naiveCenter = (minBoxPt + maxBoxPt) * 0.5;
+    T naiveRadius = 0;
+
+    for (int i = 0, icount = points.size(); i < icount; i++) {
+      const Coordinate3D<T> &point = points[i];
+
+      // Find the furthest point from the naive center to calculate the naive radius.
+      T r = (point - naiveCenter).magnitude();
+      if (r > naiveRadius) naiveRadius = r;
+
+      // Make adjustments to the Ritter Sphere to include all points.
+      T oldCenterToPointSquared = (point - ritterCenter).magnitudeSquared();
+
+      if (oldCenterToPointSquared > radiusSquared) {
+        T oldCenterToPoint = std::sqrt(oldCenterToPointSquared);
+        ritterRadius = (ritterRadius + oldCenterToPoint) * 0.5;
+
+        // Calculate center of new Ritter sphere
+        T oldToNew = oldCenterToPoint - ritterRadius;
+        ritterCenter.x = (ritterRadius * ritterCenter.x + oldToNew * point.x) / oldCenterToPoint;
+        ritterCenter.y = (ritterRadius * ritterCenter.y + oldToNew * point.y) / oldCenterToPoint;
+        ritterCenter.z = (ritterRadius * ritterCenter.z + oldToNew * point.z) / oldCenterToPoint;
+      }
+    }
+
+    // Keep the naive sphere if smaller
+    if (naiveRadius < ritterRadius) {
+      center = ritterCenter;
+      radius = ritterRadius;
+    }
+    else {
+      center = naiveCenter;
+      radius = naiveRadius;
+    }
+  }
+};
+
+/// A bounding box which is defined by a pair of minimum and maximum coordinates
+template <class T>
+class ctb::BoundingBox {
+public:
+  Coordinate3D<T> min; ///< The min coordinate of the BoundingBox
+  Coordinate3D<T> max; ///< The max coordinate of the BoundingBox
+
+  /// Create an empty BoundingBox
+  BoundingBox() {
+  }
+  /// Create a BoundingBox from the specified point stream
+  BoundingBox(const std::vector<Coordinate3D<T>> &points) {
+    fromPoints(points);
+  }
+
+  /// Calculate the BBOX from the specified point stream
+  void fromPoints(const std::vector<Coordinate3D<T>> &points) {
+    const T MAX =  std::numeric_limits<T>::infinity();
+    const T MIN = -std::numeric_limits<T>::infinity();
+    min.x = MAX;
+    min.y = MAX;
+    min.z = MAX;
+    max.x = MIN;
+    max.y = MIN;
+    max.z = MIN;
+
+    for (int i = 0, icount = points.size(); i < icount; i++) {
+      const Coordinate3D<T> &point = points[i];
+
+      if (point.x < min.x) min.x = point.x;
+      if (point.y < min.y) min.y = point.y;
+      if (point.z < min.z) min.z = point.z;
+      if (point.x > max.x) max.x = point.x;
+      if (point.y > max.y) max.y = point.y;
+      if (point.z > max.z) max.z = point.z;
+    }
+  }
+};
+
+#endif /* BBSPHERE_HPP */
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 8fca685..30be95c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -20,8 +20,14 @@ include_directories(${ZLIB_INCLUDE_DIRS})
 add_library(ctb SHARED
   GDALTile.cpp
   GDALTiler.cpp
+  GDALDatasetReader.cpp
+  CTBFileTileSerializer.cpp
+  CTBFileOutputStream.cpp
+  CTBZOutputStream.cpp
   TerrainTiler.cpp
   TerrainTile.cpp
+  MeshTiler.cpp
+  MeshTile.cpp
   GlobalMercator.cpp
   GlobalGeodetic.cpp)
 target_link_libraries(ctb ${GDAL_LIBRARIES} ${ZLIB_LIBRARIES})
@@ -29,17 +35,32 @@ target_link_libraries(ctb ${GDAL_LIBRARIES} ${ZLIB_LIBRARIES})
 # Install libctb
 set(HEADERS
   Bounds.hpp
+  BoundingSphere.hpp
   Coordinate.hpp
+  Coordinate3D.hpp
+  GDALSerializer.hpp
   GDALTile.hpp
   GDALTiler.hpp
+  GDALDatasetReader.hpp
+  CTBFileTileSerializer.hpp
+  CTBFileOutputStream.hpp
+  CTBOutputStream.hpp
+  CTBZOutputStream.hpp
   GlobalGeodetic.hpp
   GlobalMercator.hpp
   Grid.hpp
   GridIterator.hpp
+  HeightFieldChunker.hpp
+  Mesh.hpp
+  MeshIterator.hpp
+  MeshSerializer.hpp
+  MeshTile.hpp
+  MeshTiler.hpp
   RasterIterator.hpp
   RasterTiler.hpp
   CTBException.hpp
   TerrainIterator.hpp
+  TerrainSerializer.hpp
   TerrainTile.hpp
   TerrainTiler.hpp
   Tile.hpp
diff --git a/src/CTBFileOutputStream.cpp b/src/CTBFileOutputStream.cpp
new file mode 100644
index 0000000..a0dc4a3
--- /dev/null
+++ b/src/CTBFileOutputStream.cpp
@@ -0,0 +1,43 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBFileOutputStream.cpp
+ * @brief This defines the `CTBFileOutputStream` class
+ */
+
+#include "CTBFileOutputStream.hpp"
+
+using namespace ctb;
+
+/**
+ * @details 
+ * Writes a sequence of memory pointed by ptr into the FILE*.
+ */
+uint32_t
+ctb::CTBFileOutputStream::write(const void *ptr, uint32_t size) {
+  return (uint32_t)fwrite(ptr, size, 1, fp);
+}
+
+/**
+ * @details 
+ * Writes a sequence of memory pointed by ptr into the ostream. 
+ */
+uint32_t
+ctb::CTBStdOutputStream::write(const void *ptr, uint32_t size) {
+  mstream.write((const char *)ptr, size);
+  return size;
+}
diff --git a/src/CTBFileOutputStream.hpp b/src/CTBFileOutputStream.hpp
new file mode 100644
index 0000000..a1bc3bd
--- /dev/null
+++ b/src/CTBFileOutputStream.hpp
@@ -0,0 +1,60 @@
+#ifndef CTBFILEOUTPUTSTREAM_HPP
+#define CTBFILEOUTPUTSTREAM_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBFileOutputStream.hpp
+ * @brief This declares and defines the `CTBFileOutputStream` and `CTBStdOutputStream` classes
+ */
+
+#include <stdio.h>
+#include <ostream>
+#include "CTBOutputStream.hpp"
+
+namespace ctb {
+  class CTBFileOutputStream;
+  class CTBStdOutputStream;
+}
+
+/// Implements CTBOutputStream for `FILE*` objects
+class CTB_DLL ctb::CTBFileOutputStream : public ctb::CTBOutputStream {
+public:
+  CTBFileOutputStream(FILE *fptr): fp(fptr) {}
+
+  /// Writes a sequence of memory pointed by ptr into the stream
+  virtual uint32_t write(const void *ptr, uint32_t size);
+
+protected:
+  /// The underlying FILE*
+  FILE *fp;
+};
+
+/// Implements CTBOutputStream for `std::ostream` objects
+class CTB_DLL ctb::CTBStdOutputStream : public ctb::CTBOutputStream {
+public:
+  CTBStdOutputStream(std::ostream &stream): mstream(stream) {}
+
+  /// Writes a sequence of memory pointed by ptr into the stream
+  virtual uint32_t write(const void *ptr, uint32_t size);
+
+protected:
+  /// The underlying std::ostream
+  std::ostream &mstream;
+};
+
+#endif /* CTBFILEOUTPUTSTREAM_HPP */
diff --git a/src/CTBFileTileSerializer.cpp b/src/CTBFileTileSerializer.cpp
new file mode 100644
index 0000000..b067e6a
--- /dev/null
+++ b/src/CTBFileTileSerializer.cpp
@@ -0,0 +1,169 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBFileTileSerializer.cpp
+ * @brief This defines the `CTBFileTileSerializer` class
+ */
+
+#include <stdio.h>
+#include <string.h>
+#include <mutex>
+
+#include "../deps/concat.hpp"
+#include "cpl_vsi.h"
+#include "CTBException.hpp"
+#include "CTBFileTileSerializer.hpp"
+
+#include "CTBFileOutputStream.hpp"
+#include "CTBZOutputStream.hpp"
+
+using namespace std;
+using namespace ctb;
+
+#ifdef _WIN32
+static const char *osDirSep = "\\";
+#else
+static const char *osDirSep = "/";
+#endif
+
+
+/// Create a filename for a tile coordinate
+std::string
+ctb::CTBFileTileSerializer::getTileFilename(const TileCoordinate *coord, const string dirname, const char *extension) {
+  static mutex mutex;
+  VSIStatBufL stat;
+  string filename = concat(dirname, coord->zoom, osDirSep, coord->x);
+
+  lock_guard<std::mutex> lock(mutex);
+
+  // Check whether the `{zoom}/{x}` directory exists or not
+  if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) {
+    filename = concat(dirname, coord->zoom);
+
+    // Check whether the `{zoom}` directory exists or not
+    if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) {
+      // Create the `{zoom}` directory
+      if (VSIMkdir(filename.c_str(), 0755))
+        throw CTBException("Could not create the zoom level directory");
+
+    } else if (!VSI_ISDIR(stat.st_mode)) {
+      throw CTBException("Zoom level file path is not a directory");
+    }
+
+    // Create the `{zoom}/{x}` directory
+    filename += concat(osDirSep, coord->x);
+    if (VSIMkdir(filename.c_str(), 0755))
+      throw CTBException("Could not create the x level directory");
+
+  } else if (!VSI_ISDIR(stat.st_mode)) {
+    throw CTBException("X level file path is not a directory");
+  }
+
+  // Create the filename itself, adding the extension if required
+  filename += concat(osDirSep, coord->y);
+  if (extension != NULL) {
+    filename += ".";
+    filename += extension;
+  }
+
+  return filename;
+}
+
+/// Check if file exists
+static bool
+fileExists(const std::string& filename) {
+  VSIStatBufL statbuf;
+  return VSIStatExL(filename.c_str(), &statbuf, VSI_STAT_EXISTS_FLAG) == 0;
+}
+
+
+/**
+ * @details 
+ * Returns if the specified Tile Coordinate should be serialized
+ */
+bool ctb::CTBFileTileSerializer::mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) {
+  if (!mresume)
+    return true;
+
+  const string filename = getTileFilename(coordinate, moutputDir, "terrain");
+  return !fileExists(filename);
+}
+
+/**
+ * @details 
+ * Serialize a GDALTile to the Directory store
+ */
+bool 
+ctb::CTBFileTileSerializer::serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions) {
+  const TileCoordinate *coordinate = tile;
+  const string filename = getTileFilename(coordinate, moutputDir, extension);
+  const string temp_filename = concat(filename, ".tmp");
+
+  GDALDataset *poDstDS;
+  poDstDS = driver->CreateCopy(temp_filename.c_str(), tile->dataset, FALSE, creationOptions, NULL, NULL);
+
+  // Close the datasets, flushing data to destination
+  if (poDstDS == NULL) {
+    throw CTBException("Could not create GDAL tile");
+  }
+  GDALClose(poDstDS);
+
+  if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) {
+    throw CTBException("Could not rename temporary file");
+  }
+  return true;
+}
+
+/**
+ * @details 
+ * Serialize a TerrainTile to the Directory store
+ */
+bool
+ctb::CTBFileTileSerializer::serializeTile(const ctb::TerrainTile *tile) {
+  const TileCoordinate *coordinate = tile;
+  const string filename = getTileFilename(tile, moutputDir, "terrain");
+  const string temp_filename = concat(filename, ".tmp");
+
+  CTBZFileOutputStream ostream(temp_filename.c_str());
+  tile->writeFile(ostream);
+  ostream.close();
+
+  if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) {
+    throw CTBException("Could not rename temporary file");
+  }
+  return true;
+}
+
+/**
+ * @details 
+ * Serialize a MeshTile to the Directory store
+ */
+bool
+ctb::CTBFileTileSerializer::serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals) {
+  const TileCoordinate *coordinate = tile;
+  const string filename = getTileFilename(coordinate, moutputDir, "terrain");
+  const string temp_filename = concat(filename, ".tmp");
+
+  CTBZFileOutputStream ostream(temp_filename.c_str());
+  tile->writeFile(ostream, writeVertexNormals);
+  ostream.close();
+
+  if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) {
+    throw CTBException("Could not rename temporary file");
+  }
+  return true;
+}
diff --git a/src/CTBFileTileSerializer.hpp b/src/CTBFileTileSerializer.hpp
new file mode 100644
index 0000000..b723428
--- /dev/null
+++ b/src/CTBFileTileSerializer.hpp
@@ -0,0 +1,74 @@
+#ifndef CTBFILETILESERIALIZER_HPP
+#define CTBFILETILESERIALIZER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBFileTileSerializer.hpp
+ * @brief This declares and defines the `CTBFileTileSerializer` class
+ */
+
+#include <string>
+
+#include "TileCoordinate.hpp"
+#include "GDALSerializer.hpp"
+#include "TerrainSerializer.hpp"
+#include "MeshSerializer.hpp"
+
+namespace ctb {
+  class CTBFileTileSerializer;
+}
+
+/// Implements a serializer of `Tile`s based in a directory of files
+class CTB_DLL ctb::CTBFileTileSerializer : 
+  public ctb::GDALSerializer,
+  public ctb::TerrainSerializer, 
+  public ctb::MeshSerializer {
+public:
+  CTBFileTileSerializer(const std::string &outputDir, bool resume):
+    moutputDir(outputDir), 
+    mresume(resume) {}
+
+  /// Start a new serialization task
+  virtual void startSerialization() {};
+
+  /// Returns if the specified Tile Coordinate should be serialized
+  virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate);
+
+  /// Serialize a GDALTile to the store
+  virtual bool serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions);
+  /// Serialize a TerrainTile to the store
+  virtual bool serializeTile(const ctb::TerrainTile *tile);
+  /// Serialize a MeshTile to the store
+  virtual bool serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals = false);
+
+  /// Serialization finished, releases any resources loaded
+  virtual void endSerialization() {};
+
+
+  /// Create a filename for a tile coordinate
+  static std::string
+  getTileFilename(const TileCoordinate *coord, const std::string dirname, const char *extension);
+
+protected:
+  /// The target directory where serializing
+  std::string moutputDir;
+  /// Do not overwrite existing files
+  bool mresume;
+};
+
+#endif /* CTBFILETILESERIALIZER_HPP */
diff --git a/src/CTBOutputStream.hpp b/src/CTBOutputStream.hpp
new file mode 100644
index 0000000..7d0d10f
--- /dev/null
+++ b/src/CTBOutputStream.hpp
@@ -0,0 +1,40 @@
+#ifndef CTBOUTPUTSTREAM_HPP
+#define CTBOUTPUTSTREAM_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBOutputStream.hpp
+ * @brief This declares and defines the `CTBOutputStream` class
+ */
+
+#include "config.hpp"
+#include "types.hpp"
+
+namespace ctb {
+  class CTBOutputStream;
+}
+
+/// This represents a generic CTB output stream to write raw data
+class CTB_DLL ctb::CTBOutputStream {
+public:
+
+  /// Writes a sequence of memory pointed by ptr into the stream
+  virtual uint32_t write(const void *ptr, uint32_t size) = 0;
+};
+
+#endif /* CTBOUTPUTSTREAM_HPP */
diff --git a/src/CTBZOutputStream.cpp b/src/CTBZOutputStream.cpp
new file mode 100644
index 0000000..5420b5a
--- /dev/null
+++ b/src/CTBZOutputStream.cpp
@@ -0,0 +1,72 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBZOutputStream.cpp
+ * @brief This defines the `CTBZOutputStream` and `CTBZFileOutputStream` classes
+ */
+
+#include "CTBException.hpp"
+#include "CTBZOutputStream.hpp"
+
+using namespace ctb;
+
+/**
+ * @details 
+ * Writes a sequence of memory pointed by ptr into the GZFILE*.
+ */
+uint32_t
+ctb::CTBZOutputStream::write(const void *ptr, uint32_t size) {
+  if (size == 1) {
+    int c = *((const char *)ptr);
+    return gzputc(fp, c) == -1 ? 0 : 1;
+  }
+  else {
+    return gzwrite(fp, ptr, size) == 0 ? 0 : size;
+  }
+}
+
+ctb::CTBZFileOutputStream::CTBZFileOutputStream(const char *fileName) : CTBZOutputStream(NULL) {
+  gzFile file = gzopen(fileName, "wb");
+
+  if (file == NULL) {
+    throw CTBException("Failed to open file");
+  }
+  fp = file;
+}
+
+ctb::CTBZFileOutputStream::~CTBZFileOutputStream() {
+  close();
+}
+
+void
+ctb::CTBZFileOutputStream::close() {
+
+  // Try and close the file
+  if (fp) {
+    switch (gzclose(fp)) {
+    case Z_OK:
+      break;
+    case Z_STREAM_ERROR:
+    case Z_ERRNO:
+    case Z_MEM_ERROR:
+    case Z_BUF_ERROR:
+    default:
+      throw CTBException("Failed to close file");
+    }
+    fp = NULL;
+  }
+}
diff --git a/src/CTBZOutputStream.hpp b/src/CTBZOutputStream.hpp
new file mode 100644
index 0000000..b4890ba
--- /dev/null
+++ b/src/CTBZOutputStream.hpp
@@ -0,0 +1,55 @@
+#ifndef CTBZOUTPUTSTREAM_HPP
+#define CTBZOUTPUTSTREAM_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file CTBZOutputStream.hpp
+ * @brief This declares and defines the `CTBZOutputStream` class
+ */
+
+#include "zlib.h"
+#include "CTBOutputStream.hpp"
+
+namespace ctb {
+  class CTBZFileOutputStream;
+  class CTBZOutputStream;
+}
+
+/// Implements CTBOutputStream for `GZFILE` object
+class CTB_DLL ctb::CTBZOutputStream : public ctb::CTBOutputStream {
+public:
+  CTBZOutputStream(gzFile gzptr): fp(gzptr) {}
+
+  /// Writes a sequence of memory pointed by ptr into the stream
+  virtual uint32_t write(const void *ptr, uint32_t size);
+
+protected:
+  /// The underlying GZFILE*
+  gzFile fp;
+};
+
+/// Implements CTBOutputStream for gzipped files
+class CTB_DLL ctb::CTBZFileOutputStream : public ctb::CTBZOutputStream {
+public:
+  CTBZFileOutputStream(const char *fileName);
+ ~CTBZFileOutputStream();
+
+  void close();
+};
+
+#endif /* CTBZOUTPUTSTREAM_HPP */
diff --git a/src/Coordinate3D.hpp b/src/Coordinate3D.hpp
new file mode 100644
index 0000000..d4f6579
--- /dev/null
+++ b/src/Coordinate3D.hpp
@@ -0,0 +1,154 @@
+#ifndef COORDINATE3D_HPP
+#define COORDINATE3D_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+#include <vector>
+#include <cmath>
+
+/**
+ * @file Coordinate3D.hpp
+ * @brief This declares and defines the `Coordinate3D` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+namespace ctb {
+  template <class T> class Coordinate3D;
+}
+
+/// A representation of a 3-dimensional point coordinate
+template <class T>
+class ctb::Coordinate3D {
+public:
+  T x, y, z; ///< The x, y and z coordinate members
+
+  /// Create an empty coordinate
+  Coordinate3D():
+    x(0),
+    y(0),
+    z(0)
+  {}
+
+  /// The const copy constructor
+  Coordinate3D(const Coordinate3D &other):
+    x(other.x),
+    y(other.y),
+    z(other.z)
+  {}
+
+  /// Instantiate a coordinate from an x, y and z value
+  Coordinate3D(T x, T y, T z):
+    x(x),
+    y(y),
+    z(z)
+  {}
+
+  /// Overload the equality operator
+  virtual bool
+  operator==(const Coordinate3D &other) const {
+    return x == other.x
+        && y == other.y
+        && z == other.z;
+  }
+
+  /// Overload the assignment operator
+  virtual void
+  operator=(const Coordinate3D &other) {
+    x = other.x;
+    y = other.y;
+    z = other.z;
+  }
+
+  /// Gets a read-only index-ordinate of the coordinate
+  inline virtual T operator[](const int index) const {
+    return (index == 0) ? x : (index == 1 ? y : z);
+  }
+
+  /// Add operator
+  inline virtual Coordinate3D operator+(const Coordinate3D& other) const {
+    return Coordinate3D(x + other.x, y + other.y, z + other.z);
+  }
+  /// Subtract operator
+  inline virtual Coordinate3D operator-(const Coordinate3D& other) const {
+    return Coordinate3D(x - other.x, y - other.y, z - other.z);
+  }
+  /// Multiply operator
+  inline virtual Coordinate3D operator*(const Coordinate3D& other) const {
+    return Coordinate3D(x * other.x, y * other.y, z * other.z);
+  }
+  /// Divide operator
+  inline virtual Coordinate3D operator/(const Coordinate3D& other) const {
+    return Coordinate3D(x / other.x, y / other.y, z / other.z);
+  }
+
+  /// AddByScalar operator
+  inline virtual Coordinate3D operator+(const T scalar) const {
+    return Coordinate3D(x + scalar, y + scalar, z + scalar);
+  }
+  /// SubtractByScalar operator
+  inline virtual Coordinate3D operator-(const T scalar) const {
+    return Coordinate3D(x - scalar, y - scalar, z - scalar);
+  }
+  /// MultiplyByScalar operator
+  inline virtual Coordinate3D operator*(const T scalar) const {
+    return Coordinate3D(x * scalar, y * scalar, z * scalar);
+  }
+  /// DivideByScalar operator
+  inline virtual Coordinate3D operator/(const T scalar) const {
+    return Coordinate3D(x / scalar, y / scalar, z / scalar);
+  }
+
+  /// Cross product
+  inline Coordinate3D<T> cross(const Coordinate3D<T> &other) const {
+    return Coordinate3D((y * other.z) - (other.y * z), 
+                        (z * other.x) - (other.z * x), 
+                        (x * other.y) - (other.x * y));
+  }
+  /// Dot product
+  inline double dot(const Coordinate3D<T> &other) const {
+    return (x * other.x) + (y * other.y) + (z * other.z);
+  }
+
+  // Cartesian3d methods
+  inline T magnitudeSquared(void) const {
+    return (x * x) + (y * y) + (z * z);
+  }
+  inline T magnitude(void) const {
+    return std::sqrt(magnitudeSquared());
+  }
+  inline static Coordinate3D<T> add(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    return p1 + p2;
+  }
+  inline static Coordinate3D<T> subtract(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    return p1 - p2;
+  }
+  inline static T distanceSquared(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    T xdiff = p1.x - p2.x;
+    T ydiff = p1.y - p2.y;
+    T zdiff = p1.z - p2.z;
+    return (xdiff * xdiff) + (ydiff * ydiff) + (zdiff * zdiff);
+  }
+  inline static T distance(const Coordinate3D<T> &p1, const Coordinate3D<T> &p2) {
+    return std::sqrt(distanceSquared(p1, p2));
+  }
+  inline Coordinate3D<T> normalize(void) const {
+    T mgn = magnitude();
+    return Coordinate3D(x / mgn, y / mgn, z / mgn);
+  }
+};
+
+#endif /* COORDINATE3D_HPP */
diff --git a/src/GDALDatasetReader.cpp b/src/GDALDatasetReader.cpp
new file mode 100644
index 0000000..43f97ef
--- /dev/null
+++ b/src/GDALDatasetReader.cpp
@@ -0,0 +1,156 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file GDALDatasetReader.cpp
+ * @brief This defines the `GDALDatasetReader` class
+ */
+
+#include "gdal_priv.h"
+#include "gdalwarper.h"
+
+#include "CTBException.hpp"
+#include "GDALDatasetReader.hpp"
+#include "TerrainTiler.hpp"
+
+using namespace ctb;
+
+/**
+ * @details 
+ * Read a region of raster heights for the specified Dataset and Coordinate. 
+ * This method uses `GDALRasterBand::RasterIO` function.
+ */
+float *
+ctb::GDALDatasetReader::readRasterHeights(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) {
+  GDALTile *rasterTile = createRasterTile(tiler, dataset, coord); // the raster associated with this tile coordinate
+
+  const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY;
+  float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float));
+
+  GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1);
+
+  if (heightsBand->RasterIO(GF_Read, 0, 0, tileSizeX, tileSizeY,
+                            (void *) rasterHeights, tileSizeX, tileSizeY, GDT_Float32,
+                            0, 0) != CE_None) {
+    delete rasterTile;
+    CPLFree(rasterHeights);
+
+    throw CTBException("Could not read heights from raster");
+  }
+  delete rasterTile;
+  return rasterHeights;
+}
+
+/// Create a raster tile from a tile coordinate
+GDALTile *
+ctb::GDALDatasetReader::createRasterTile(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord) {
+  return tiler.createRasterTile(dataset, coord);
+}
+
+/// Create a VTR raster overview from a GDALDataset
+GDALDataset *
+ctb::GDALDatasetReader::createOverview(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, int overviewIndex) {
+  int nFactorScale = 2 << overviewIndex;
+  int nRasterXSize = dataset->GetRasterXSize() / nFactorScale;
+  int nRasterYSize = dataset->GetRasterYSize() / nFactorScale;
+
+  GDALDataset *poOverview = NULL;
+  double adfGeoTransform[6];
+
+  // Should We create an overview of the Dataset?
+  if (nRasterXSize > 4 && nRasterYSize > 4 && dataset->GetGeoTransform(adfGeoTransform) == CE_None) {
+    adfGeoTransform[1] *= nFactorScale;
+    adfGeoTransform[5] *= nFactorScale;
+
+    TerrainTiler tempTiler(tiler.dataset(), tiler.grid(), tiler.options);
+    tempTiler.crsWKT = "";
+    GDALTile *rasterTile = createRasterTile(tempTiler, dataset, coord);
+    if (rasterTile) {
+      poOverview = rasterTile->detach();
+      delete rasterTile;
+    }
+  }
+  return poOverview;
+}
+
+/// The destructor
+ctb::GDALDatasetReaderWithOverviews::~GDALDatasetReaderWithOverviews() {
+  reset();
+}
+
+/// Read a region of raster heights into an array for the specified Dataset and Coordinate
+float *
+ctb::GDALDatasetReaderWithOverviews::readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) {
+  GDALDataset *mainDataset = dataset;
+
+  const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY;
+  float *rasterHeights = (float *)CPLCalloc(TILE_CELL_SIZE, sizeof(float));
+
+  // Replace GDAL Dataset by last valid Overview.
+  for (int i = mOverviews.size() - 1; i >= 0; --i) {
+    if (mOverviews[i]) {
+      dataset = mOverviews[i];
+      break;
+    }
+  }
+
+  // Extract the raster data, using overviews when necessary
+  bool rasterOk = false;
+  while (!rasterOk) {
+    GDALTile *rasterTile = createRasterTile(poTiler, dataset, coord); // the raster associated with this tile coordinate
+
+    GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1);
+
+    if (heightsBand->RasterIO(GF_Read, 0, 0, tileSizeX, tileSizeY,
+                              (void *) rasterHeights, tileSizeX, tileSizeY, GDT_Float32,
+                              0, 0) != CE_None) {
+      
+      GDALDataset *psOverview = createOverview(poTiler, mainDataset, coord, mOverviewIndex++);
+      if (psOverview) {
+        mOverviews.push_back(psOverview);
+        dataset = psOverview;
+      }
+      else {
+        delete rasterTile;
+        CPLFree(rasterHeights);
+        throw CTBException("Could not create an overview of current GDAL dataset");
+      }
+    }
+    else {
+      rasterOk = true;
+    }
+    delete rasterTile;
+  }
+
+  // Everything ok?
+  if (!rasterOk) {
+    CPLFree(rasterHeights);
+    throw CTBException("Could not read heights from raster");
+  }
+  return rasterHeights;
+}
+
+/// Releases all overviews
+void 
+ctb::GDALDatasetReaderWithOverviews::reset() {
+  mOverviewIndex = 0;
+
+  for (int i = mOverviews.size() - 1; i >= 0; --i) {
+    GDALDataset *poOverview = mOverviews[i];
+    GDALClose(poOverview);
+  }
+  mOverviews.clear();
+}
diff --git a/src/GDALDatasetReader.hpp b/src/GDALDatasetReader.hpp
new file mode 100644
index 0000000..a3ef18c
--- /dev/null
+++ b/src/GDALDatasetReader.hpp
@@ -0,0 +1,100 @@
+#ifndef GDALDATASETREADER_HPP
+#define GDALDATASETREADER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file GDALDatasetReader.hpp
+ * @brief This declares the `GDALDatasetReader` class
+ */
+
+#include <string>
+#include <vector>
+#include "gdalwarper.h"
+
+#include "TileCoordinate.hpp"
+#include "GDALTiler.hpp"
+
+namespace ctb {
+  class GDALDatasetReader;
+  class GDALDatasetReaderWithOverviews;
+}
+
+/**
+ * @brief Read raster tiles from a GDAL Dataset
+ *
+ * This abstract base class is associated with a GDAL dataset. 
+ * It allows to read a region of the raster according to
+ * a region defined by a Tile Coordinate.
+ *
+ * We can define our own
+ */
+class CTB_DLL ctb::GDALDatasetReader {
+public:
+  /// Read a region of raster heights into an array for the specified Dataset and Coordinate
+  static float *
+  readRasterHeights(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY);
+
+  /// Read a region of raster heights into an array for the specified Dataset and Coordinate
+  virtual float *
+  readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) = 0;
+
+protected:
+  /// Create a raster tile from a tile coordinate
+  static GDALTile *
+  createRasterTile(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord);
+
+  /// Create a VTR raster overview from a GDALDataset
+  static GDALDataset *
+  createOverview(const GDALTiler &tiler, GDALDataset *dataset, const TileCoordinate &coord, int overviewIndex);
+};
+
+/** 
+ * @brief Implements a GDALDatasetReader that takes care of 'Integer overflow' errors.
+ * 
+ * This class creates Overviews to avoid 'Integer overflow' errors when extracting 
+ * raster data.
+ */
+class CTB_DLL ctb::GDALDatasetReaderWithOverviews : public ctb::GDALDatasetReader {
+public:
+
+  /// Instantiate a GDALDatasetReaderWithOverviews
+  GDALDatasetReaderWithOverviews(const GDALTiler &tiler): 
+    poTiler(tiler), 
+    mOverviewIndex(0) {}
+
+  /// The destructor
+  ~GDALDatasetReaderWithOverviews();
+
+  /// Read a region of raster heights into an array for the specified Dataset and Coordinate
+  virtual float *
+  readRasterHeights(GDALDataset *dataset, const TileCoordinate &coord, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) override;
+
+  /// Releases all overviews
+  void reset();
+
+protected:
+  /// The tiler to use
+  const GDALTiler &poTiler;
+
+  /// List of VRT Overviews of the underlying GDAL dataset
+  std::vector<GDALDataset *> mOverviews;
+  /// Current VRT Overview
+  int mOverviewIndex;
+};
+
+#endif /* GDALDATASETREADER_HPP */
diff --git a/src/GDALSerializer.hpp b/src/GDALSerializer.hpp
new file mode 100644
index 0000000..130d8c5
--- /dev/null
+++ b/src/GDALSerializer.hpp
@@ -0,0 +1,52 @@
+#ifndef GDALSERIALIZER_HPP
+#define GDALSERIALIZER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file GDALSerializer.hpp
+ * @brief This declares and defines the `GDALSerializer` class
+ */
+
+#include "config.hpp"
+#include "TileCoordinate.hpp"
+#include "GDALTile.hpp"
+
+#include "cpl_string.h"
+
+namespace ctb {
+  class GDALSerializer;
+}
+
+/// Store `GDALTile`s from a GDAL Dataset
+class CTB_DLL ctb::GDALSerializer {
+public:
+
+  /// Start a new serialization task
+  virtual void startSerialization() = 0;
+
+  /// Returns if the specified Tile Coordinate should be serialized
+  virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0;
+
+  /// Serialize a GDALTile to the store
+  virtual bool serializeTile(const ctb::GDALTile *tile, GDALDriver *driver, const char *extension, CPLStringList &creationOptions) = 0;
+
+  /// Serialization finished, releases any resources loaded
+  virtual void endSerialization() = 0;
+};
+
+#endif /* GDALSERIALIZER_HPP */
diff --git a/src/GDALTile.cpp b/src/GDALTile.cpp
index ad7df74..c3a6dae 100644
--- a/src/GDALTile.cpp
+++ b/src/GDALTile.cpp
@@ -34,3 +34,18 @@ GDALTile::~GDALTile() {
     }
   }
 }
+
+/// Detach the underlying GDAL dataset
+GDALDataset *GDALTile::detach() {
+  if (dataset != NULL) {
+    GDALDataset *poDataset = dataset;
+    dataset = NULL;
+
+    if (transformer != NULL) {
+      GDALDestroyGenImgProjTransformer(transformer);
+      transformer = NULL;
+    }
+    return poDataset;
+  }
+  return NULL;
+}
diff --git a/src/GDALTile.hpp b/src/GDALTile.hpp
index 86ba270..85914b7 100644
--- a/src/GDALTile.hpp
+++ b/src/GDALTile.hpp
@@ -57,6 +57,9 @@ class CTB_DLL ctb::GDALTile :
 
   GDALDataset *dataset;
 
+  /// Detach the underlying GDAL dataset
+  GDALDataset *detach();
+
 protected:
   friend class GDALTiler;
 
diff --git a/src/GDALTiler.cpp b/src/GDALTiler.cpp
index 70fc735..0fb216e 100644
--- a/src/GDALTiler.cpp
+++ b/src/GDALTiler.cpp
@@ -70,6 +70,11 @@ GDALTiler::GDALTiler(GDALDataset *poDataset, const Grid &grid, const TilerOption
     OGRSpatialReference srcSRS = OGRSpatialReference(srcWKT);
     OGRSpatialReference gridSRS = mGrid.getSRS();
 
+    #if ( GDAL_VERSION_MAJOR >= 3 )
+    srcSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+    gridSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+    #endif
+
     if (!srcSRS.IsSame(&gridSRS)) { // it doesn't match
       // Check the srs is valid
       switch(srcSRS.Validate()) {
@@ -173,7 +178,7 @@ GDALTiler::~GDALTiler() {
 }
 
 GDALTile *
-GDALTiler::createRasterTile(const TileCoordinate &coord) const {
+GDALTiler::createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const {
   // Convert the tile bounds into a geo transform
   double adfGeoTransform[6],
     resolution = mGrid.resolution(coord.zoom);
@@ -186,7 +191,7 @@ GDALTiler::createRasterTile(const TileCoordinate &coord) const {
   adfGeoTransform[4] = 0;
   adfGeoTransform[5] = -resolution;
 
-  GDALTile *tile = createRasterTile(adfGeoTransform);
+  GDALTile *tile = createRasterTile(dataset, adfGeoTransform);
   static_cast<TileCoordinate &>(*tile) = coord;
 
   // Set the shifted geo transform to the VRT
@@ -209,6 +214,12 @@ GDALTiler::createRasterTile(const TileCoordinate &coord) const {
  * This code is adapted from that found in `gdalwarp.cpp` implementing the
  * `gdalwarp -ovr` option.
  */
+#if ( GDAL_VERSION_MAJOR >= 3 )
+#include "gdaloverviewdataset.cpp"
+#elif ( GDAL_VERSION_MAJOR >= 2 && GDAL_VERSION_MINOR >= 2 )
+#include "gdaloverviewdataset-gdal2x.cpp"
+#endif
+
 static
 GDALDatasetH
 getOverviewDataset(GDALDatasetH hSrcDS, GDALTransformerFunc pfnTransformer, void *hTransformerArg) {
@@ -246,7 +257,11 @@ getOverviewDataset(GDALDatasetH hSrcDS, GDALTransformerFunc pfnTransformer, void
               if( iOvr >= 0 )
                 {
                   //std::cout << "CTB WARPING: Selecting overview level " << iOvr << " for output dataset " << nPixels << "x" << nLines << std::endl;
-                  poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE);
+                #if ( GDAL_VERSION_MAJOR >= 3 || ( GDAL_VERSION_MAJOR >= 2 && GDAL_VERSION_MINOR >= 2 ) )
+                  poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE );
+                #else
+                  poSrcOvrDS = GDALCreateOverviewDataset( poSrcDS, iOvr, FALSE, FALSE );
+                #endif
                 }
             }
         }
@@ -267,13 +282,13 @@ getOverviewDataset(GDALDatasetH hSrcDS, GDALTransformerFunc pfnTransformer, void
  * dataset.
  */
 GDALTile *
-GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const {
-  if (poDataset == NULL) {
+GDALTiler::createRasterTile(GDALDataset *dataset, double (&adfGeoTransform)[6]) const {
+  if (dataset == NULL) {
     throw CTBException("No GDAL dataset is set");
   }
 
   // The source and sink datasets
-  GDALDatasetH hSrcDS = (GDALDatasetH) dataset();
+  GDALDatasetH hSrcDS = (GDALDatasetH) dataset;
   GDALDatasetH hDstDS;
 
   // The transformation option list
@@ -304,7 +319,25 @@ GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const {
   psWarpOptions->panDstBands =
     (int *) CPLMalloc(sizeof(int) * psWarpOptions->nBandCount );
 
+  psWarpOptions->padfSrcNoDataReal =
+    (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));
+  psWarpOptions->padfSrcNoDataImag =
+    (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));
+  psWarpOptions->padfDstNoDataReal =
+    (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));
+  psWarpOptions->padfDstNoDataImag =
+    (double *)CPLCalloc(psWarpOptions->nBandCount, sizeof(double));
+
   for (short unsigned int i = 0; i < psWarpOptions->nBandCount; ++i) {
+    int bGotNoData = FALSE;
+    double noDataValue = poDataset->GetRasterBand(i + 1)->GetNoDataValue(&bGotNoData);
+    if (!bGotNoData) noDataValue = -32768;
+
+    psWarpOptions->padfSrcNoDataReal[i] = noDataValue;
+    psWarpOptions->padfSrcNoDataImag[i] = 0;
+    psWarpOptions->padfDstNoDataReal[i] = noDataValue;
+    psWarpOptions->padfDstNoDataImag[i] = 0;
+
     psWarpOptions->panDstBands[i] = psWarpOptions->panSrcBands[i] = i + 1;
   }
 
@@ -315,23 +348,29 @@ GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const {
     throw CTBException("Could not create image to image transformer");
   }
 
+  // Specify the destination geotransform
+  GDALSetGenImgProjTransformerDstGeoTransform(transformerArg, adfGeoTransform );
+
   // Try and get an overview from the source dataset that corresponds more
   // closely to the resolution of this tile.
   GDALDatasetH hWrkSrcDS = getOverviewDataset(hSrcDS, GDALGenImgProjTransform, transformerArg);
   if (hWrkSrcDS == NULL) {
     hWrkSrcDS = psWarpOptions->hSrcDS = hSrcDS;
   } else {
+    psWarpOptions->hSrcDS = hWrkSrcDS;
+
     // We need to recreate the transform when operating on an overview.
     GDALDestroyGenImgProjTransformer( transformerArg );
+
     transformerArg = GDALCreateGenImgProjTransformer2( hWrkSrcDS, NULL, transformOptions.List() );
     if(transformerArg == NULL) {
       GDALDestroyWarpOptions(psWarpOptions);
       throw CTBException("Could not create overview image to image transformer");
     }
-  }
 
-  // Specify the destination geotransform
-  GDALSetGenImgProjTransformerDstGeoTransform(transformerArg, adfGeoTransform );
+    // Specify the destination geotransform
+    GDALSetGenImgProjTransformerDstGeoTransform(transformerArg, adfGeoTransform );
+  }
 
   // Decide if we are doing an approximate or exact transformation
   if (options.errorThreshold) {
@@ -353,11 +392,6 @@ GDALTiler::createRasterTile(double (&adfGeoTransform)[6]) const {
     psWarpOptions->pfnTransformer = GDALGenImgProjTransform;
   }
 
-  // Specify a multi threaded warp operation using all CPU cores
-  CPLStringList warpOptions(psWarpOptions->papszWarpOptions, false);
-  warpOptions.SetNameValue("NUM_THREADS", "ALL_CPUS");
-  psWarpOptions->papszWarpOptions = warpOptions.StealList();
-
   // The raster tile is represented as a VRT dataset
   hDstDS = GDALCreateWarpedVRT(hWrkSrcDS, mGrid.tileSize(), mGrid.tileSize(), adfGeoTransform, psWarpOptions);
 
diff --git a/src/GDALTiler.hpp b/src/GDALTiler.hpp
index 1f3a3ee..87477db 100644
--- a/src/GDALTiler.hpp
+++ b/src/GDALTiler.hpp
@@ -33,6 +33,7 @@
 namespace ctb {
   struct TilerOptions;
   class GDALTiler;
+  class GDALDatasetReader; // forward declaration
 }
 
 /// Options passed to a `GDALTiler`
@@ -91,7 +92,7 @@ class CTB_DLL ctb::GDALTiler {
 
   /// Create a tile from a tile coordinate
   virtual Tile *
-  createTile(const TileCoordinate &coord) const = 0;
+  createTile(GDALDataset *dataset, const TileCoordinate &coord) const = 0;
 
   /// Get the maximum zoom level for the dataset
   inline i_zoom
@@ -151,16 +152,18 @@ class CTB_DLL ctb::GDALTiler {
   }
 
 protected:
+  friend class GDALDatasetReader;
+
   /// Close the underlying dataset
   void closeDataset();
 
   /// Create a raster tile from a tile coordinate
   virtual GDALTile *
-  createRasterTile(const TileCoordinate &coord) const;
+  createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const;
 
   /// Create a raster tile from a geo transform
   virtual GDALTile *
-  createRasterTile(double (&adfGeoTransform)[6]) const;
+  createRasterTile(GDALDataset *dataset, double (&adfGeoTransform)[6]) const;
 
   /// The grid used for generating tiles
   Grid mGrid;
diff --git a/src/GlobalGeodetic.cpp b/src/GlobalGeodetic.cpp
index 5a34016..c741d2b 100644
--- a/src/GlobalGeodetic.cpp
+++ b/src/GlobalGeodetic.cpp
@@ -27,6 +27,11 @@ using namespace ctb;
 static OGRSpatialReference
 setSRS(void) {
   OGRSpatialReference srs;
+  
+  #if ( GDAL_VERSION_MAJOR >= 3 )
+  srs.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+  #endif
+  
   srs.importFromEPSG(4326);
   return srs;
 }
diff --git a/src/GlobalMercator.cpp b/src/GlobalMercator.cpp
index d808bb7..9c50b1f 100644
--- a/src/GlobalMercator.cpp
+++ b/src/GlobalMercator.cpp
@@ -34,6 +34,11 @@ const double GlobalMercator::cOriginShift = GlobalMercator::cEarthCircumference
 static OGRSpatialReference
 setSRS(void) {
   OGRSpatialReference srs;
+  
+  #if ( GDAL_VERSION_MAJOR >= 3 )
+  srs.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+  #endif
+  
   srs.importFromEPSG(3857);
   return srs;
 }
diff --git a/src/Grid.hpp b/src/Grid.hpp
index 156e23e..9702370 100644
--- a/src/Grid.hpp
+++ b/src/Grid.hpp
@@ -67,7 +67,11 @@ class ctb::Grid {
     mXOriginShift(extent.getWidth() / 2),
     mYOriginShift(extent.getHeight() / 2),
     mZoomFactor(zoomFactor)
-  {}
+  {
+    #if ( GDAL_VERSION_MAJOR >= 3 )
+    mSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+    #endif
+  }
 
   /// Overload the assignment operator
   Grid &
diff --git a/src/HeightFieldChunker.hpp b/src/HeightFieldChunker.hpp
new file mode 100644
index 0000000..9bab8ff
--- /dev/null
+++ b/src/HeightFieldChunker.hpp
@@ -0,0 +1,493 @@
+#ifndef HEIGHTFIELDCHUNKER_HPP
+#define HEIGHTFIELDCHUNKER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file HeightFieldChucker.hpp
+ * @brief This declares and defines the `mesh` and `heightfield` classes
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include <vector>
+
+#include "cpl_config.h"
+#include "cpl_string.h"
+
+/**
+ * Helper classes to fill an irregular mesh of triangles from a heightmap tile.
+ * They are a refactored version from 'heightfield_chunker.cpp' from 
+ * http://tulrich.com/geekstuff/chunklod.html
+ *
+ * It applies the Chunked LOD strategy by 'Thatcher Ulrich'
+ * preserving the input geometric error.
+ */
+namespace ctb { namespace chunk {
+  struct gen_state;
+  class mesh;
+  class heightfield;
+} }
+
+/// Helper struct with state info for chunking a HeightField.
+struct ctb::chunk::gen_state {
+  int my_buffer[2][2];  // x,y coords of the last two vertices emitted by the generate_ functions.
+  int activation_level; // for determining whether a vertex is enabled in the block we're working on
+  int ptr;              // indexes my_buffer.
+  int previous_level;   // for keeping track of level changes during recursion.
+
+  /// Returns true if the specified vertex is in my_buffer.
+  bool in_my_buffer(int x, int y) const
+  {
+    return ((x == my_buffer[0][0]) && (y == my_buffer[0][1]))
+        || ((x == my_buffer[1][0]) && (y == my_buffer[1][1]));
+  }
+
+  /// Sets the current my_buffer entry to (x,y)
+  void set_my_buffer(int x, int y)
+  {
+    my_buffer[ptr][0] = x;
+    my_buffer[ptr][1] = y;
+  }
+};
+
+/// An irregular mesh of triangles target of the HeightField chunker process.
+class ctb::chunk::mesh {
+public:
+  /// Clear all data.
+  virtual void clear() = 0;
+
+  /// New vertex (Call this in strip order).
+  virtual void emit_vertex(const heightfield &heightfield, int x, int y) = 0;
+};
+
+/// Defines a regular grid of heigths or HeightField.
+class ctb::chunk::heightfield {
+public:
+  /// Constructor
+  heightfield(float *tileHeights, int tileSize) {
+    int tileCellSize = tileSize * tileSize;
+
+    m_heights = tileHeights;
+    m_size = tileSize;
+    m_log_size = (int)(log2((float)m_size - 1) + 0.5);
+
+    // Initialize level array.
+    m_levels = (int*)CPLMalloc(tileCellSize * sizeof(int));
+    for (int i = 0; i < tileCellSize; i++) m_levels[i] = 255;
+  }
+  ~heightfield() {
+    clear();
+  }
+
+  /// Apply the specified maximum geometric error to fill the level info of the grid.
+  void applyGeometricError(double maximumGeometricError, bool smoothSmallZooms = false) {
+    int tileCellSize = m_size * m_size;
+
+    // Initialize level array.
+    for (int i = 0; i < tileCellSize; i++) m_levels[i] = 255;
+
+    // Run a view-independent L-K style BTT update on the heightfield,
+    // to generate error and activation_level values for each element.
+    update(maximumGeometricError, 0, m_size - 1, m_size - 1, m_size - 1, 0, 0); // sw half of the square
+    update(maximumGeometricError, m_size - 1, 0, 0, 0, m_size - 1, m_size - 1); // ne half of the square
+
+    // Make sure our corner verts are activated.
+    int size = (m_size - 1);
+    activate(size, 0, 0);
+    activate(0, 0, 0);
+    activate(0, size, 0);
+    activate(size, size, 0);
+
+    // Activate some vertices to smooth the shape of the Globe for small zooms.
+    if (smoothSmallZooms) {
+      int step = size / 16;
+
+      for (int x = 0; x <= size; x += step) {
+        for (int y = 0; y <= size; y += step) {
+          if (get_level(x, y) == -1) activate(x, y, 0);
+        }
+      }
+    }
+
+    // Propagate the activation_level values of verts to their parent verts,
+    // quadtree LOD style. Gives same result as L-K.
+    for (int i = 0; i < m_log_size; i++) {
+      propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i);
+      propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i);
+    }
+  }
+
+  /// Returns the Coordinate of the Neighbor of the specified border (Left=0, Top=1, Right=2, Botton=3).
+  static ctb::TileCoordinate neighborCoord(const Grid& grid, const ctb::TileCoordinate &coord, int borderIndex, bool& okNeighborCoord) {
+    okNeighborCoord = true;
+
+    switch (borderIndex)
+    {
+    case 0:
+      if (coord.x <= 0) {
+        okNeighborCoord = false;
+        return TileCoordinate();
+      }
+      return ctb::TileCoordinate(coord.zoom, coord.x - 1, coord.y);
+
+    case 1:
+      if (coord.y >= grid.getTileExtent(coord.zoom).getMaxY()) {
+        okNeighborCoord = false;
+        return TileCoordinate();
+      }
+      return ctb::TileCoordinate(coord.zoom, coord.x, coord.y + 1);
+
+    case 2: 
+      if (coord.x >= grid.getTileExtent(coord.zoom).getMaxX()) {
+        okNeighborCoord = false;
+        return TileCoordinate();
+      }
+      return ctb::TileCoordinate(coord.zoom, coord.x + 1, coord.y);
+
+    case 3: 
+      if (coord.y <= 0) {
+        okNeighborCoord = false;
+        return TileCoordinate();
+      }
+      return ctb::TileCoordinate(coord.zoom, coord.x, coord.y - 1);
+
+    default:
+      throw CTBException("Bad Neighbor border index");
+    }
+  }
+
+  /// Apply the activation state of the border of the specified Neighbor.
+  void applyBorderActivationState(const ctb::chunk::heightfield &hf, int borderIndex) {
+    int level = -1;
+
+    switch (borderIndex) //-> (Left=0, Top=1, Right=2, Botton=3)
+    {
+    case 0:
+      for (int x = m_size - 1, y = 0; y < m_size; y++) {
+        level = hf.get_level(x, y);
+        if (level != -1) activate(0, y, level);
+      }
+      break;
+    case 1:
+      for (int x = 0, y = m_size - 1; x < m_size; x++) {
+        level = hf.get_level(x, y);
+        if (level != -1) activate(x, 0, level);
+      }
+      break;
+    case 2:
+      for (int x = 0, y = 0; y < m_size; y++) {
+        level = hf.get_level(x, y);
+        if (level != -1) activate(m_size - 1, y, level);
+      }
+      break;
+    case 3:
+      for (int x = 0, y = 0; x < m_size; x++) {
+        level = hf.get_level(x, y);
+        if (level != -1) activate(x, m_size - 1, level);
+      }
+      break;
+    default:
+      throw CTBException("Bad Neighbor border index");
+    }
+
+    // Propagate the activation_level values of verts to their parent verts,
+    // quadtree LOD style. Gives same result as L-K.
+    for (int i = 0; i < m_log_size; i++) {
+      propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i);
+      propagate_activation_level(m_size >> 1, m_size >> 1, m_log_size - 1, i);
+    }
+  }
+
+  /// Clear all object data
+  void clear() {
+    m_heights = NULL;
+    m_size = 0;
+    m_log_size = 0;
+
+    if (m_levels) {
+      CPLFree(m_levels);
+      m_levels = NULL;
+    }
+  }
+
+  /// Return the array-index of specified coordinate, row order by default.
+  virtual int indexOfGridCoordinate(int x, int y) const {
+    return (y * m_size) + x;
+  }
+  /// Return the height of specified coordinate.
+  virtual float height(int x, int y) const {
+    int index = indexOfGridCoordinate(x, y);
+    return m_heights[index];
+  }
+
+  /// Generates the mesh using verts which are active at the given level.
+  void generateMesh(ctb::chunk::mesh &mesh, int level) {
+    int x0 = 0;
+    int y0 = 0;
+
+    int size = (1 << m_log_size);
+    int half_size = size >> 1;
+    int cx = x0 + half_size;
+    int cy = y0 + half_size;
+
+    // Start making the mesh.
+    mesh.clear();
+
+    // !!! This needs to be done in propagate, or something (too late now) !!!
+    // Make sure our corner verts are activated on this level.
+    activate(x0 + size, y0, level);
+    activate(x0, y0, level);
+    activate(x0, y0 + size, level);
+    activate(x0 + size, y0 + size, level);
+
+    // Generate the mesh.
+    const heightfield &hf = *this;
+    generate_block(hf, mesh, level, m_log_size, x0 + half_size, y0 + half_size);
+  }
+
+private:
+  int m_size;         // Number of cols and rows of this Heightmap
+  int m_log_size;     // size == (1 << log_size) + 1
+  float *m_heights;   // grid of heights
+  int *m_levels;      // grid of activation levels
+
+  /// Return the activation level at (x, y)
+  int get_level(int x, int y) const
+  {
+    int index = indexOfGridCoordinate(x, y);
+    int level = m_levels[index];
+
+    if (x & 1) {
+      level = level >> 4;
+    }
+    level &= 0x0F;
+    if (level == 0x0F) return -1;
+    else return level;
+  }
+  /// Set the activation level at (x, y)
+  void set_level(int x, int y, int newlevel)
+  {
+    newlevel &= 0x0F;
+    int index = indexOfGridCoordinate(x, y);
+    int level = m_levels[index];
+
+    if (x & 1) {
+      level = (level & 0x0F) | (newlevel << 4);
+    }
+    else {
+      level = (level & 0xF0) | (newlevel);
+    }
+    m_levels[index] = level;
+  }
+  /// Sets the activation_level to the given level.
+  /// if it's greater than the vert's current activation level.
+  void activate(int x, int y, int level)
+  {
+    int current_level = get_level(x, y);
+    if (level > current_level) set_level(x, y, level);
+  }
+
+  /// Given the triangle, computes an error value and activation level
+  /// for its base vertex, and recurses to child triangles.
+  bool update(double base_max_error, int ax, int ay, int rx, int ry, int lx, int ly)
+  {
+    bool res = false;
+
+    // Compute the coordinates of this triangle's base vertex.
+    int dx = lx - rx;
+    int dy = ly - ry;
+
+    if (std::abs(dx) <= 1 && std::abs(dy) <= 1) {
+      // We've reached the base level.  There's no base
+      // vertex to update, and no child triangles to
+      // recurse to.
+
+      return false;
+    }
+
+    // base vert is midway between left and right verts.
+    int bx = rx + (dx >> 1);
+    int by = ry + (dy >> 1);
+
+    float heightB = height(bx, by);
+    float heightL = height(lx, ly);
+    float heightR = height(rx, ry);
+    float error_B = std::abs(heightB - 0.5 * (heightL + heightR));
+
+    if (error_B >= base_max_error) {
+      // Compute the mesh level above which this vertex
+      // needs to be included in LOD meshes.
+      int activation_level = (int)std::floor(log2(error_B / base_max_error) + 0.5);
+
+      // Force the base vert to at least this activation level.
+      activate(bx, by, activation_level);
+      res = true;
+    }
+
+    // Recurse to child triangles.
+    update(base_max_error, bx, by, ax, ay, rx, ry); // base, apex, right
+    update(base_max_error, bx, by, lx, ly, ax, ay); // base, left, apex
+
+    return res;
+  }
+
+  /// Does a quadtree descent through the heightfield, in the square with
+  /// center at (cx, cz) and size of (2 ^ (level + 1) + 1).  Descends
+  /// until level == target_level, and then propagates this square's
+  /// child center verts to the corresponding edge vert, and the edge
+  /// verts to the center.  Essentially the quadtree meshing update
+  /// dependency graph as in my Gamasutra article.  Must call this with
+  /// successively increasing target_level to get correct propagation.
+  void propagate_activation_level(int cx, int cy, int level, int target_level)
+  {
+    int half_size = 1 << level;
+    int quarter_size = half_size >> 1;
+
+    if (level > target_level) {
+      // Recurse to children.
+      for (int j = 0; j < 2; j++) {
+        for (int i = 0; i < 2; i++) {
+          propagate_activation_level(
+            cx - quarter_size + half_size * i,
+            cy - quarter_size + half_size * j,
+            level - 1, target_level);
+        }
+      }
+      return;
+    }
+
+    // We're at the target level. Do the propagation on this square.
+    if (level > 0) {
+      int lev = 0;
+
+      // Propagate child verts to edge verts.
+      lev = get_level(cx + quarter_size, cy - quarter_size); // ne.
+      activate(cx + half_size, cy, lev);
+      activate(cx, cy - half_size, lev);
+
+      lev = get_level(cx - quarter_size, cy - quarter_size); // nw.
+      activate(cx, cy - half_size, lev);
+      activate(cx - half_size, cy, lev);
+
+      lev = get_level(cx - quarter_size, cy + quarter_size); // sw.
+      activate(cx - half_size, cy, lev);
+      activate(cx, cy + half_size, lev);
+
+      lev = get_level(cx + quarter_size, cy + quarter_size); // se.
+      activate(cx, cy + half_size, lev);
+      activate(cx + half_size, cy, lev);
+    }
+
+    // Propagate edge verts to center.
+    activate(cx, cy, get_level(cx + half_size, cy));
+    activate(cx, cy, get_level(cx, cy - half_size));
+    activate(cx, cy, get_level(cx, cy + half_size));
+    activate(cx, cy, get_level(cx - half_size, cy));
+  }
+
+  /// Auxiliary function for generate_block().
+  /// Generates a mesh from a triangular quadrant of a square heightfield block.
+  /// Paraphrased directly out of Lindstrom et al, SIGGRAPH '96.
+  void generate_quadrant(const heightfield &hf, mesh &mesh, gen_state* state, int lx, int ly, int tx, int ty, int rx, int ry, int recursion_level) const {
+    if (recursion_level <= 0) return;
+
+    if (hf.get_level(tx, ty) >= state->activation_level) {
+      // Find base vertex.
+      int bx = (lx + rx) >> 1;
+      int by = (ly + ry) >> 1;
+
+      generate_quadrant(hf, mesh, state, lx, ly, bx, by, tx, ty, recursion_level - 1); // left half of quadrant
+
+      if (state->in_my_buffer(tx, ty) == false) {
+        if ((recursion_level + state->previous_level) & 1) {
+          state->ptr ^= 1;
+        }
+        else {
+          int x = state->my_buffer[1 - state->ptr][0];
+          int y = state->my_buffer[1 - state->ptr][1];
+          mesh.emit_vertex(hf, x, y); // or, emit vertex(last - 1);
+        }
+        mesh.emit_vertex(hf, tx, ty);
+        state->set_my_buffer(tx, ty);
+        state->previous_level = recursion_level;
+      }
+      generate_quadrant(hf, mesh, state, tx, ty, bx, by, rx, ry, recursion_level - 1);
+    }
+  }
+  /// Generate the mesh for the specified square with the given center.
+  /// This is paraphrased directly out of Lindstrom et al, SIGGRAPH '96.
+  /// It generates a square mesh by walking counterclockwise around four
+  /// triangular quadrants.
+  /// The resulting mesh is composed of a single continuous triangle strip,
+  /// with a few corners turned via degenerate tris where necessary.
+  void generate_block(const heightfield &hf, mesh &mesh, int activation_level, int log_size, int cx, int cy) const {
+    int hs = 1 << (log_size - 1);
+
+    // quadrant corner coordinates.
+    int q[4][2] = {
+      { cx + hs, cy + hs }, // se
+      { cx + hs, cy - hs }, // ne
+      { cx - hs, cy - hs }, // nw
+      { cx - hs, cy + hs }, // sw
+    };
+
+    // Init state for generating mesh.
+    gen_state state;
+    state.ptr = 0;
+    state.previous_level = 0;
+    state.activation_level = activation_level;
+    for (int i = 0; i < 4; i++) {
+      state.my_buffer[i >> 1][i & 1] = -1;
+    }
+
+    mesh.emit_vertex(hf,q[0][0], q[0][1]);
+    state.set_my_buffer(q[0][0], q[0][1]);
+
+    {for (int i = 0; i < 4; i++) {
+      if ((state.previous_level & 1) == 0) {
+        // tulrich: turn a corner?
+        state.ptr ^= 1;
+      }
+      else {
+        // tulrich: jump via degenerate?
+        int x = state.my_buffer[1 - state.ptr][0];
+        int y = state.my_buffer[1 - state.ptr][1];
+
+        mesh.emit_vertex(hf, x, y); // or, emit vertex(last - 1);
+      }
+
+      // Initial vertex of quadrant.
+      mesh.emit_vertex(hf,q[i][0], q[i][1]);
+      state.set_my_buffer(q[i][0], q[i][1]);
+      state.previous_level = 2 * log_size + 1;
+
+      generate_quadrant(hf, mesh,
+        &state,
+        q[i][0], q[i][1], // q[i][l]
+        cx, cy, // q[i][t]
+        q[(i + 1) & 3][0], q[(i + 1) & 3][1], // q[i][r]
+        2 * log_size
+      );
+    }}
+    if (state.in_my_buffer(q[0][0], q[0][1]) == false) {
+      // finish off the strip.  @@ may not be necessary?
+      mesh.emit_vertex(hf, q[0][0], q[0][1]);
+    }
+  }
+};
+
+#endif /* HEIGHTFIELDCHUNKER_HPP */
diff --git a/src/Mesh.hpp b/src/Mesh.hpp
new file mode 100644
index 0000000..2e3559a
--- /dev/null
+++ b/src/Mesh.hpp
@@ -0,0 +1,82 @@
+#ifndef CTBMESH_HPP
+#define CTBMESH_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file Mesh.hpp
+ * @brief This declares the `Mesh` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include <cstdint>
+#include <vector>
+#include <cstring>
+#include "types.hpp"
+#include "CTBException.hpp"
+
+namespace ctb {
+  class Mesh;
+}
+
+/**
+ * @brief An abstract base class for a mesh of triangles
+ */
+class ctb::Mesh
+{
+public:
+
+  /// Create an empty mesh
+  Mesh()
+  {}
+
+public:
+
+  /// The array of shared vertices of a mesh
+  std::vector<CRSVertex> vertices;
+
+  /// The index collection for each triangle in the mesh (3 for each triangle)
+  std::vector<uint32_t> indices;
+
+  /// Write mesh data to a WKT file
+  void writeWktFile(const char *fileName) const {
+    FILE *fp = fopen(fileName, "w");
+
+    if (fp == NULL) {
+      throw CTBException("Failed to open file");
+    }
+
+    char wktText[512];
+    memset(wktText, 0, sizeof(wktText));
+
+    for (int i = 0, icount = indices.size(); i < icount; i += 3) {
+      CRSVertex v0 = vertices[indices[i]];
+      CRSVertex v1 = vertices[indices[i+1]];
+      CRSVertex v2 = vertices[indices[i+2]];
+
+      sprintf(wktText, "(%.8f %.8f %f, %.8f %.8f %f, %.8f %.8f %f, %.8f %.8f %f)",
+        v0.x, v0.y, v0.z,
+        v1.x, v1.y, v1.z,
+        v2.x, v2.y, v2.z,
+        v0.x, v0.y, v0.z);
+      fprintf(fp, "POLYGON Z(%s)\n", wktText);
+    }
+    fclose(fp);
+  };
+};
+
+#endif /* CTBMESH_HPP */
diff --git a/src/MeshIterator.hpp b/src/MeshIterator.hpp
new file mode 100644
index 0000000..3383b67
--- /dev/null
+++ b/src/MeshIterator.hpp
@@ -0,0 +1,72 @@
+#ifndef MESHITERATOR_HPP
+#define MESHITERATOR_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshIterator.hpp
+ * @brief This declares the `MeshIterator` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include "MeshTiler.hpp"
+#include "GridIterator.hpp"
+
+namespace ctb {
+  class MeshIterator;
+}
+
+/**
+ * @brief This forward iterates over all `MeshTile`s in a `MeshTiler`
+ *
+ * Instances of this class take a `MeshTiler` in the constructor and are used
+ * to forward iterate over all tiles in the tiler, returning a `MeshTile *`
+ * when dereferenced.  It is the caller's responsibility to call `delete` on the
+ * tile.
+ */
+class ctb::MeshIterator :
+  public GridIterator
+{
+public:
+
+  /// Instantiate an iterator with a tiler
+  MeshIterator(const MeshTiler &tiler) :
+    MeshIterator(tiler, tiler.maxZoomLevel(), 0)
+  {}
+
+  MeshIterator(const MeshTiler &tiler, i_zoom startZoom, i_zoom endZoom = 0) :
+    GridIterator(tiler.grid(), tiler.bounds(), startZoom, endZoom),
+    tiler(tiler)
+  {}
+
+  /// Override the dereference operator to return a Tile
+  virtual MeshTile *
+  operator*() const {
+    return tiler.createMesh(tiler.dataset(), *(GridIterator::operator*()));
+  }
+
+  virtual MeshTile *
+  operator*(ctb::GDALDatasetReader *reader) const {
+    return tiler.createMesh(tiler.dataset(), *(GridIterator::operator*()), reader);
+  }
+
+protected:
+
+  const MeshTiler &tiler;              ///< The tiler we are iterating over
+};
+
+#endif /* MESHITERATOR_HPP */
diff --git a/src/MeshSerializer.hpp b/src/MeshSerializer.hpp
new file mode 100644
index 0000000..98d6fa0
--- /dev/null
+++ b/src/MeshSerializer.hpp
@@ -0,0 +1,50 @@
+#ifndef MESHSERIALIZER_HPP
+#define MESHSERIALIZER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshSerializer.hpp
+ * @brief This declares and defines the `MeshSerializer` class
+ */
+
+#include "config.hpp"
+#include "TileCoordinate.hpp"
+#include "MeshTile.hpp"
+
+namespace ctb {
+  class MeshSerializer;
+}
+
+/// Store `MeshTile`s from a GDAL Dataset
+class CTB_DLL ctb::MeshSerializer {
+public:
+
+  /// Start a new serialization task
+  virtual void startSerialization() = 0;
+
+  /// Returns if the specified Tile Coordinate should be serialized
+  virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0;
+
+  /// Serialize a MeshTile to the store
+  virtual bool serializeTile(const ctb::MeshTile *tile, bool writeVertexNormals = false) = 0;
+
+  /// Serialization finished, releases any resources loaded
+  virtual void endSerialization() = 0;
+};
+
+#endif /* MESHSERIALIZER_HPP */
diff --git a/src/MeshTile.cpp b/src/MeshTile.cpp
new file mode 100644
index 0000000..7173dfb
--- /dev/null
+++ b/src/MeshTile.cpp
@@ -0,0 +1,447 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshTile.cpp
+ * @brief This defines the `MeshTile` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include <cmath>
+#include <vector>
+#include <map>
+#include "cpl_conv.h"
+
+#include "CTBException.hpp"
+#include "MeshTile.hpp"
+#include "BoundingSphere.hpp"
+#include "CTBZOutputStream.hpp"
+
+using namespace ctb;
+
+////////////////////////////////////////////////////////////////////////////////
+// Utility functions
+
+// Constants taken from http://cesiumjs.org/2013/04/25/Horizon-culling
+double llh_ecef_radiusX = 6378137.0;
+double llh_ecef_radiusY = 6378137.0;
+double llh_ecef_radiusZ = 6356752.3142451793;
+
+double llh_ecef_rX = 1.0 / llh_ecef_radiusX;
+double llh_ecef_rY = 1.0 / llh_ecef_radiusY;
+double llh_ecef_rZ = 1.0 / llh_ecef_radiusZ;
+
+// Stolen from https://github.com/bistromath/gr-air-modes/blob/master/python/mlat.py
+// WGS84 reference ellipsoid constants
+// http://en.wikipedia.org/wiki/Geodetic_datum#Conversion_calculations
+// http://en.wikipedia.org/wiki/File%3aECEF.png
+//
+double llh_ecef_wgs84_a = llh_ecef_radiusX;       // Semi - major axis
+double llh_ecef_wgs84_b = llh_ecef_radiusZ;       // Semi - minor axis
+double llh_ecef_wgs84_e2 = 0.0066943799901975848; // First eccentricity squared
+
+// LLH2ECEF
+static inline double llh_ecef_n(double x) {
+  double snx = std::sin(x);
+  return llh_ecef_wgs84_a / std::sqrt(1.0 - llh_ecef_wgs84_e2 * (snx * snx));
+}
+static inline CRSVertex LLH2ECEF(const CRSVertex& coordinate) {
+  double lon = coordinate.x * (M_PI / 180.0);
+  double lat = coordinate.y * (M_PI / 180.0);
+  double alt = coordinate.z;
+
+  double x = (llh_ecef_n(lat) + alt) * std::cos(lat) * std::cos(lon);
+  double y = (llh_ecef_n(lat) + alt) * std::cos(lat) * std::sin(lon);
+  double z = (llh_ecef_n(lat) * (1.0 - llh_ecef_wgs84_e2) + alt) * std::sin(lat);
+
+  return CRSVertex(x, y, z);
+}
+
+// HORIZON OCCLUSION POINT
+// https://cesiumjs.org/2013/05/09/Computing-the-horizon-occlusion-point
+//
+static inline double ocp_computeMagnitude(const CRSVertex &position, const CRSVertex &sphereCenter) {
+  double magnitudeSquared = position.magnitudeSquared();
+  double magnitude = std::sqrt(magnitudeSquared);
+  CRSVertex direction = position * (1.0 / magnitude);
+
+  // For the purpose of this computation, points below the ellipsoid
+  // are considered to be on it instead.
+  magnitudeSquared = std::fmax(1.0, magnitudeSquared);
+  magnitude = std::fmax(1.0, magnitude);
+
+  double cosAlpha = direction.dot(sphereCenter);
+  double sinAlpha = direction.cross(sphereCenter).magnitude();
+  double cosBeta = 1.0 / magnitude;
+  double sinBeta = std::sqrt(magnitudeSquared - 1.0) * cosBeta;
+
+  return 1.0 / (cosAlpha * cosBeta - sinAlpha * sinBeta);
+}
+static inline CRSVertex ocp_fromPoints(const std::vector<CRSVertex> &points, const BoundingSphere<double> &boundingSphere) {
+  const double MIN = -std::numeric_limits<double>::infinity();
+  double max_magnitude = MIN;
+
+  // Bring coordinates to ellipsoid scaled coordinates
+  const CRSVertex &center = boundingSphere.center;
+  CRSVertex scaledCenter = CRSVertex(center.x * llh_ecef_rX, center.y * llh_ecef_rY, center.z * llh_ecef_rZ);
+
+  for (int i = 0, icount = points.size(); i < icount; i++) {
+    const CRSVertex &point = points[i];
+    CRSVertex scaledPoint(point.x * llh_ecef_rX, point.y * llh_ecef_rY, point.z * llh_ecef_rZ);
+
+    double magnitude = ocp_computeMagnitude(scaledPoint, scaledCenter);
+    if (magnitude > max_magnitude) max_magnitude = magnitude;
+  }
+  return scaledCenter * max_magnitude;
+}
+
+// PACKAGE IO
+const double SHORT_MAX = 32767.0;
+const int BYTESPLIT = 65636;
+
+static inline int quantizeIndices(const double &origin, const double &factor, const double &value) {
+  return int(std::round((value - origin) * factor));
+}
+
+// Write the edge indices of the mesh
+template <typename T> int writeEdgeIndices(CTBOutputStream &ostream, const Mesh &mesh, double edgeCoord, int componentIndex) {
+  std::vector<uint32_t> indices;
+  std::map<uint32_t, size_t> ihash;
+
+  for (size_t i = 0, icount = mesh.indices.size(); i < icount; i++) {
+    uint32_t indice = mesh.indices[i];
+    double val = mesh.vertices[indice][componentIndex];
+
+    if (val == edgeCoord) {
+      std::map<uint32_t, size_t>::iterator it = ihash.find(indice);
+
+      if (it == ihash.end()) {
+        ihash.insert(std::make_pair(indice, i));
+        indices.push_back(indice);
+      }
+    }
+  }
+
+  int edgeCount = indices.size();
+  ostream.write(&edgeCount, sizeof(int));
+
+  for (size_t i = 0; i < edgeCount; i++) {
+    T indice = (T)indices[i];
+    ostream.write(&indice, sizeof(T));
+  }
+  return indices.size();
+}
+
+// ZigZag-Encodes a number (-1 = 1, -2 = 3, 0 = 0, 1 = 2, 2 = 4)
+static inline uint16_t zigZagEncode(int n) {
+  return (n << 1) ^ (n >> 31);
+}
+
+// Triangle area
+static inline double triangleArea(const CRSVertex &a, const CRSVertex &b) {
+  double i = std::pow(a[1] * b[2] - a[2] * b[1], 2);
+  double j = std::pow(a[2] * b[0] - a[0] * b[2], 2);
+  double k = std::pow(a[0] * b[1] - a[1] * b[0], 2);
+  return 0.5 * sqrt(i + j + k);
+}
+
+// Constraint a value to lie between two values
+static inline double clamp_value(double value, double min, double max) {
+  return value < min ? min : value > max ? max : value;
+}
+
+// Converts a scalar value in the range [-1.0, 1.0] to a SNORM in the range [0, rangeMax]
+static inline unsigned char snorm_value(double value, double rangeMax = 255) {
+  return (unsigned char)int(std::round((clamp_value(value, -1.0, 1.0) * 0.5 + 0.5) * rangeMax));
+}
+
+/**
+ * Encodes a normalized vector into 2 SNORM values in the range of [0-rangeMax] following the 'oct' encoding.
+ *
+ * Oct encoding is a compact representation of unit length vectors.
+ * The 'oct' encoding is described in "A Survey of Efficient Representations of Independent Unit Vectors",
+ * Cigolle et al 2014: {@link http://jcgt.org/published/0003/02/01/}
+ */
+static inline Coordinate<unsigned char> octEncode(const CRSVertex &vector, double rangeMax = 255) {
+  Coordinate<double> temp;
+  double llnorm = std::abs(vector.x) + std::abs(vector.y) + std::abs(vector.z);
+  temp.x = vector.x / llnorm;
+  temp.y = vector.y / llnorm;
+
+  if (vector.z < 0) {
+    double x = temp.x;
+    double y = temp.y;
+    temp.x = (1.0 - std::abs(y)) * (x < 0.0 ? -1.0 : 1.0);
+    temp.y = (1.0 - std::abs(x)) * (y < 0.0 ? -1.0 : 1.0);
+  }
+  return Coordinate<unsigned char>(snorm_value(temp.x, rangeMax), snorm_value(temp.y, rangeMax));
+}
+
+////////////////////////////////////////////////////////////////////////////////
+
+MeshTile::MeshTile():
+  Tile(),
+  mChildren(0)
+{}
+
+MeshTile::MeshTile(const TileCoordinate &coord):
+  Tile(coord),
+  mChildren(0)
+{}
+
+/**
+ * @details This writes gzipped terrain data to a file.
+ */
+void 
+MeshTile::writeFile(const char *fileName, bool writeVertexNormals) const {
+  CTBZFileOutputStream ostream(fileName);
+  writeFile(ostream);
+}
+
+/**
+ * @details This writes raw terrain data to an output stream.
+ */
+void
+MeshTile::writeFile(CTBOutputStream &ostream, bool writeVertexNormals) const {
+
+  // Calculate main header mesh data
+  std::vector<CRSVertex> cartesianVertices;
+  BoundingSphere<double> cartesianBoundingSphere;
+  BoundingBox<double> cartesianBounds;
+  BoundingBox<double> bounds;
+
+  cartesianVertices.resize(mMesh.vertices.size());
+  for (size_t i = 0, icount = mMesh.vertices.size(); i < icount; i++) {
+    const CRSVertex &vertex = mMesh.vertices[i];
+    cartesianVertices[i] = LLH2ECEF(vertex);
+  }
+  cartesianBoundingSphere.fromPoints(cartesianVertices);
+  cartesianBounds.fromPoints(cartesianVertices);
+  bounds.fromPoints(mMesh.vertices);
+
+
+  // # Write the mesh header data:
+  // # https://github.com/AnalyticalGraphicsInc/quantized-mesh
+  //
+  // The center of the tile in Earth-centered Fixed coordinates.
+  double centerX = cartesianBounds.min.x + 0.5 * (cartesianBounds.max.x - cartesianBounds.min.x);
+  double centerY = cartesianBounds.min.y + 0.5 * (cartesianBounds.max.y - cartesianBounds.min.y);
+  double centerZ = cartesianBounds.min.z + 0.5 * (cartesianBounds.max.z - cartesianBounds.min.z);
+  ostream.write(&centerX, sizeof(double));
+  ostream.write(&centerY, sizeof(double));
+  ostream.write(&centerZ, sizeof(double));
+  //
+  // The minimum and maximum heights in the area covered by this tile.
+  float minimumHeight = (float)bounds.min.z;
+  float maximumHeight = (float)bounds.max.z;
+  ostream.write(&minimumHeight, sizeof(float));
+  ostream.write(&maximumHeight, sizeof(float));
+  //
+  // The tile's bounding sphere. The X,Y,Z coordinates are again expressed
+  // in Earth-centered Fixed coordinates, and the radius is in meters.
+  double boundingSphereCenterX = cartesianBoundingSphere.center.x;
+  double boundingSphereCenterY = cartesianBoundingSphere.center.y;
+  double boundingSphereCenterZ = cartesianBoundingSphere.center.z;
+  double boundingSphereRadius  = cartesianBoundingSphere.radius;
+  ostream.write(&boundingSphereCenterX, sizeof(double));
+  ostream.write(&boundingSphereCenterY, sizeof(double));
+  ostream.write(&boundingSphereCenterZ, sizeof(double));
+  ostream.write(&boundingSphereRadius , sizeof(double));
+  //
+  // The horizon occlusion point, expressed in the ellipsoid-scaled Earth-centered Fixed frame.
+  CRSVertex horizonOcclusionPoint = ocp_fromPoints(cartesianVertices, cartesianBoundingSphere);
+  ostream.write(&horizonOcclusionPoint.x, sizeof(double));
+  ostream.write(&horizonOcclusionPoint.y, sizeof(double));
+  ostream.write(&horizonOcclusionPoint.z, sizeof(double));
+
+
+  // # Write mesh vertices (X Y Z components of each vertex):
+  int vertexCount = mMesh.vertices.size();
+  ostream.write(&vertexCount, sizeof(int));
+  for (int c = 0; c < 3; c++) {
+    double origin = bounds.min[c];
+    double factor = 0;
+    if (bounds.max[c] > bounds.min[c]) factor = SHORT_MAX / (bounds.max[c] - bounds.min[c]);
+
+    // Move the initial value
+    int u0 = quantizeIndices(origin, factor, mMesh.vertices[0][c]), u1, ud;
+    uint16_t sval = zigZagEncode(u0);
+    ostream.write(&sval, sizeof(uint16_t));
+
+    for (size_t i = 1, icount = mMesh.vertices.size(); i < icount; i++) {
+      u1 = quantizeIndices(origin, factor, mMesh.vertices[i][c]);
+      ud = u1 - u0;
+      sval = zigZagEncode(ud);
+      ostream.write(&sval, sizeof(uint16_t));
+      u0 = u1;
+    }
+  }
+
+  // # Write mesh indices:
+  int triangleCount = mMesh.indices.size() / 3;
+  ostream.write(&triangleCount, sizeof(int));
+  if (vertexCount > BYTESPLIT) {
+    uint32_t highest = 0;
+    uint32_t code;
+
+    // Write main indices
+    for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) {
+      code = highest - mMesh.indices[i];
+      ostream.write(&code, sizeof(uint32_t));
+      if (code == 0) highest++;
+    }
+
+    // Write all vertices on the edge of the tile (W, S, E, N)
+    writeEdgeIndices<uint32_t>(ostream, mMesh, bounds.min.x, 0);
+    writeEdgeIndices<uint32_t>(ostream, mMesh, bounds.min.y, 1);
+    writeEdgeIndices<uint32_t>(ostream, mMesh, bounds.max.x, 0);
+    writeEdgeIndices<uint32_t>(ostream, mMesh, bounds.max.y, 1);
+  }
+  else {
+    uint16_t highest = 0;
+    uint16_t code;
+
+    // Write main indices
+    for (size_t i = 0, icount = mMesh.indices.size(); i < icount; i++) {
+      code = highest - mMesh.indices[i];
+      ostream.write(&code, sizeof(uint16_t));
+      if (code == 0) highest++;
+    }
+
+    // Write all vertices on the edge of the tile (W, S, E, N)
+    writeEdgeIndices<uint16_t>(ostream, mMesh, bounds.min.x, 0);
+    writeEdgeIndices<uint16_t>(ostream, mMesh, bounds.min.y, 1);
+    writeEdgeIndices<uint16_t>(ostream, mMesh, bounds.max.x, 0);
+    writeEdgeIndices<uint16_t>(ostream, mMesh, bounds.max.y, 1);
+  }
+
+  // # Write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting:
+  if (writeVertexNormals && triangleCount > 0) {
+    unsigned char extensionId = 1;
+    ostream.write(&extensionId, sizeof(unsigned char));
+    int extensionLength = 2 * vertexCount;
+    ostream.write(&extensionLength, sizeof(int));
+
+    std::vector<CRSVertex> normalsPerVertex(vertexCount);
+    std::vector<CRSVertex> normalsPerFace(triangleCount);
+    std::vector<double> areasPerFace(triangleCount);
+
+    for (size_t i = 0, icount = mMesh.indices.size(), j = 0; i < icount; i+=3, j++) {
+      const CRSVertex &v0 = cartesianVertices[ mMesh.indices[i  ] ];
+      const CRSVertex &v1 = cartesianVertices[ mMesh.indices[i+1] ];
+      const CRSVertex &v2 = cartesianVertices[ mMesh.indices[i+2] ];
+
+      CRSVertex normal = (v1 - v0).cross(v2 - v0);
+      double area = triangleArea(v0, v1);
+      normalsPerFace[j] = normal;
+      areasPerFace[j] = area;
+    }
+    for (size_t i = 0, icount = mMesh.indices.size(), j = 0; i < icount; i+=3, j++) {
+      int indexV0 = mMesh.indices[i  ];
+      int indexV1 = mMesh.indices[i+1];
+      int indexV2 = mMesh.indices[i+2];
+
+      CRSVertex weightedNormal = normalsPerFace[j] * areasPerFace[j];
+
+      normalsPerVertex[indexV0] = normalsPerVertex[indexV0] + weightedNormal;
+      normalsPerVertex[indexV1] = normalsPerVertex[indexV1] + weightedNormal;
+      normalsPerVertex[indexV2] = normalsPerVertex[indexV2] + weightedNormal;
+    }
+    for (size_t i = 0; i < vertexCount; i++) {
+      Coordinate<unsigned char> xy = octEncode(normalsPerVertex[i].normalize());
+      ostream.write(&xy.x, sizeof(unsigned char));
+      ostream.write(&xy.y, sizeof(unsigned char));
+    }
+  }
+}
+
+bool
+MeshTile::hasChildren() const {
+  return mChildren;
+}
+
+bool
+MeshTile::hasChildSW() const {
+  return ((mChildren & TERRAIN_CHILD_SW) == TERRAIN_CHILD_SW);
+}
+
+bool
+MeshTile::hasChildSE() const {
+  return ((mChildren & TERRAIN_CHILD_SE) == TERRAIN_CHILD_SE);
+}
+
+bool
+MeshTile::hasChildNW() const {
+  return ((mChildren & TERRAIN_CHILD_NW) == TERRAIN_CHILD_NW);
+}
+
+bool
+MeshTile::hasChildNE() const {
+  return ((mChildren & TERRAIN_CHILD_NE) == TERRAIN_CHILD_NE);
+}
+
+void
+MeshTile::setChildSW(bool on) {
+  if (on) {
+    mChildren |= TERRAIN_CHILD_SW;
+  } else {
+    mChildren &= ~TERRAIN_CHILD_SW;
+  }
+}
+
+void
+MeshTile::setChildSE(bool on) {
+  if (on) {
+    mChildren |= TERRAIN_CHILD_SE;
+  } else {
+    mChildren &= ~TERRAIN_CHILD_SE;
+  }
+}
+
+void
+MeshTile::setChildNW(bool on) {
+  if (on) {
+    mChildren |= TERRAIN_CHILD_NW;
+  } else {
+    mChildren &= ~TERRAIN_CHILD_NW;
+  }
+}
+
+void
+MeshTile::setChildNE(bool on) {
+  if (on) {
+    mChildren |= TERRAIN_CHILD_NE;
+  } else {
+    mChildren &= ~TERRAIN_CHILD_NE;
+  }
+}
+
+void
+MeshTile::setAllChildren(bool on) {
+  if (on) {
+    mChildren = TERRAIN_CHILD_SW | TERRAIN_CHILD_SE | TERRAIN_CHILD_NW | TERRAIN_CHILD_NE;
+  } else {
+    mChildren = 0;
+  }
+}
+
+const Mesh & MeshTile::getMesh() const {
+  return mMesh;
+}
+
+Mesh & MeshTile::getMesh() {
+  return mMesh;
+}
diff --git a/src/MeshTile.hpp b/src/MeshTile.hpp
new file mode 100644
index 0000000..2300c37
--- /dev/null
+++ b/src/MeshTile.hpp
@@ -0,0 +1,132 @@
+#ifndef MESHTILE_HPP
+#define MESHTILE_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshTile.hpp
+ * @brief This declares the `MeshTile` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include "config.hpp"
+#include "Mesh.hpp"
+#include "TileCoordinate.hpp"
+#include "Tile.hpp"
+#include "CTBOutputStream.hpp"
+
+namespace ctb {
+  class MeshTile;
+}
+
+/**
+ * @brief `Terrain` data associated with a `Mesh`
+ *
+ * This aims to implement the Cesium [quantized-mesh-1.0 terrain
+ * format](https://github.com/AnalyticalGraphicsInc/quantized-mesh).
+ */
+class CTB_DLL ctb::MeshTile :
+  public Tile
+{
+  friend class MeshTiler;
+
+public:
+
+  /// Create an empty mesh tile object
+  MeshTile();
+
+  /// Create a mesh tile from a tile coordinate
+  MeshTile(const TileCoordinate &coord);
+
+  /// Write terrain data to the filesystem
+  void
+  writeFile(const char *fileName, bool writeVertexNormals = false) const;
+
+  /// Write terrain data to an output stream
+  void
+  writeFile(CTBOutputStream &ostream, bool writeVertexNormals = false) const;
+
+  /// Does the terrain tile have child tiles?
+  bool
+  hasChildren() const;
+
+  /// Does the terrain tile have a south west child tile?
+  bool
+  hasChildSW() const;
+
+  /// Does the terrain tile have a south east child tile?
+  bool
+  hasChildSE() const;
+
+  /// Does the terrain tile have a north west child tile?
+  bool
+  hasChildNW() const;
+
+  /// Does the terrain tile have a north east child tile?
+  bool
+  hasChildNE() const;
+
+  /// Specify that there is a south west child tile
+  void
+  setChildSW(bool on = true);
+
+  /// Specify that there is a south east child tile
+  void
+  setChildSE(bool on = true);
+
+  /// Specify that there is a north west child tile
+  void
+  setChildNW(bool on = true);
+
+  /// Specify that there is a north east child tile
+  void
+  setChildNE(bool on = true);
+
+  /// Specify that all child tiles are present
+  void
+  setAllChildren(bool on = true);
+
+  /// Get the mesh data as a const object
+  const ctb::Mesh & getMesh() const;
+
+  /// Get the mesh data
+  ctb::Mesh & getMesh();
+
+protected:
+
+  /// The terrain mesh data
+  ctb::Mesh mMesh;
+
+private:
+
+  char mChildren;               ///< The child flags
+
+  /**
+   * @brief Bit flags defining child tile existence
+   *
+   * There is a good discussion on bitflags
+   * [here](http://www.dylanleigh.net/notes/c-cpp-tricks.html#Using_"Bitflags").
+   */
+  enum Children {
+    TERRAIN_CHILD_SW = 1,       // 2^0, bit 0
+    TERRAIN_CHILD_SE = 2,       // 2^1, bit 1
+    TERRAIN_CHILD_NW = 4,       // 2^2, bit 2
+    TERRAIN_CHILD_NE = 8        // 2^3, bit 3
+  };
+};
+
+#endif /* MESHTILE_HPP */
diff --git a/src/MeshTiler.cpp b/src/MeshTiler.cpp
new file mode 100644
index 0000000..88192bd
--- /dev/null
+++ b/src/MeshTiler.cpp
@@ -0,0 +1,236 @@
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshTiler.cpp
+ * @brief This defines the `MeshTiler` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include "CTBException.hpp"
+#include "MeshTiler.hpp"
+#include "HeightFieldChunker.hpp"
+#include "GDALDatasetReader.hpp"
+
+using namespace ctb;
+
+////////////////////////////////////////////////////////////////////////////////
+
+/**
+ * Implementation of ctb::chunk::mesh for ctb::Mesh class.
+ */
+class WrapperMesh : public ctb::chunk::mesh {
+private:
+  CRSBounds &mBounds;
+  Mesh &mMesh;
+  double mCellSizeX;
+  double mCellSizeY;
+
+  std::map<int, int> mIndicesMap;
+  Coordinate<int> mTriangles[3];
+  bool mTriOddOrder;
+  int mTriIndex;
+
+public:
+  WrapperMesh(CRSBounds &bounds, Mesh &mesh, i_tile tileSizeX, i_tile tileSizeY):
+    mMesh(mesh),
+    mBounds(bounds),
+    mTriOddOrder(false),
+    mTriIndex(0) {
+    mCellSizeX = (bounds.getMaxX() - bounds.getMinX()) / (double)(tileSizeX - 1);
+    mCellSizeY = (bounds.getMaxY() - bounds.getMinY()) / (double)(tileSizeY - 1);
+  }
+
+  virtual void clear() {
+    mMesh.vertices.clear();
+    mMesh.indices.clear();
+    mIndicesMap.clear();
+    mTriOddOrder = false;
+    mTriIndex = 0;
+  }
+  virtual void emit_vertex(const ctb::chunk::heightfield &heightfield, int x, int y) {
+    mTriangles[mTriIndex].x = x;
+    mTriangles[mTriIndex].y = y;
+    mTriIndex++;
+
+    if (mTriIndex == 3) {
+      mTriOddOrder = !mTriOddOrder;
+
+      if (mTriOddOrder) {
+        appendVertex(heightfield, mTriangles[0].x, mTriangles[0].y);
+        appendVertex(heightfield, mTriangles[1].x, mTriangles[1].y);
+        appendVertex(heightfield, mTriangles[2].x, mTriangles[2].y);
+      }
+      else {
+        appendVertex(heightfield, mTriangles[1].x, mTriangles[1].y);
+        appendVertex(heightfield, mTriangles[0].x, mTriangles[0].y);
+        appendVertex(heightfield, mTriangles[2].x, mTriangles[2].y);
+      }
+      mTriangles[0].x = mTriangles[1].x;
+      mTriangles[0].y = mTriangles[1].y;
+      mTriangles[1].x = mTriangles[2].x;
+      mTriangles[1].y = mTriangles[2].y;
+      mTriIndex--;
+    }
+  }
+  void appendVertex(const ctb::chunk::heightfield &heightfield, int x, int y) {
+    int iv;
+    int index = heightfield.indexOfGridCoordinate(x, y);
+
+    std::map<int, int>::iterator it = mIndicesMap.find(index);
+
+    if (it == mIndicesMap.end()) {
+      iv = mMesh.vertices.size();
+
+      double xmin = mBounds.getMinX();
+      double ymax = mBounds.getMaxY();
+      double height = heightfield.height(x, y);
+
+      mMesh.vertices.push_back(CRSVertex(xmin + (x * mCellSizeX), ymax - (y * mCellSizeY), height));
+      mIndicesMap.insert(std::make_pair(index, iv));
+    }
+    else {
+      iv = it->second;
+    }
+    mMesh.indices.push_back(iv);
+  }
+};
+
+////////////////////////////////////////////////////////////////////////////////
+
+void 
+ctb::MeshTiler::prepareSettingsOfTile(MeshTile *terrainTile, GDALDataset *dataset, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const {
+  const ctb::i_tile TILE_SIZE = tileSizeX;
+
+  // Number of tiles in the horizontal direction at tile level zero.
+  double resolutionAtLevelZero = mGrid.resolution(0);
+  int numberOfTilesAtLevelZero = (int)(mGrid.getExtent().getWidth() / (tileSizeX * resolutionAtLevelZero));
+  // Default quality of terrain created from heightmaps (TerrainProvider.js).
+  double heightmapTerrainQuality = 0.25;
+  // Earth semi-major-axis in meters.
+  const double semiMajorAxis = 6378137.0;
+  // Appropriate geometric error estimate when the geometry comes from a heightmap (TerrainProvider.js).
+  double maximumGeometricError = MeshTiler::getEstimatedLevelZeroGeometricErrorForAHeightmap(
+    semiMajorAxis,
+    heightmapTerrainQuality * mMeshQualityFactor,
+    TILE_SIZE,
+    numberOfTilesAtLevelZero
+  );
+  // Geometric error for current Level.
+  maximumGeometricError /= (double)(1 << coord.zoom);
+
+  // Convert the raster grid into an irregular mesh applying the Chunked LOD strategy by 'Thatcher Ulrich'.
+  // http://tulrich.com/geekstuff/chunklod.html
+  //
+  ctb::chunk::heightfield heightfield(rasterHeights, TILE_SIZE);
+  heightfield.applyGeometricError(maximumGeometricError, coord.zoom <= 6);
+  //
+  // Propagate the geometric error of neighbours to avoid gaps in borders.
+  if (coord.zoom > 6) {
+    ctb::CRSBounds datasetBounds = bounds();
+
+    for (int borderIndex = 0; borderIndex < 4; borderIndex++) {
+      bool okNeighborCoord = true;
+      ctb::TileCoordinate neighborCoord = ctb::chunk::heightfield::neighborCoord(mGrid, coord, borderIndex, okNeighborCoord);
+      if (!okNeighborCoord)
+        continue;
+
+      ctb::CRSBounds neighborBounds = mGrid.tileBounds(neighborCoord);
+
+      if (datasetBounds.overlaps(neighborBounds)) {
+        float *neighborHeights = ctb::GDALDatasetReader::readRasterHeights(*this, dataset, neighborCoord, mGrid.tileSize(), mGrid.tileSize());
+
+        ctb::chunk::heightfield neighborHeightfield(neighborHeights, TILE_SIZE);
+        neighborHeightfield.applyGeometricError(maximumGeometricError);
+        heightfield.applyBorderActivationState(neighborHeightfield, borderIndex);
+
+        CPLFree(neighborHeights);
+      }
+    }
+  }
+  ctb::CRSBounds mGridBounds = mGrid.tileBounds(coord);
+  Mesh &tileMesh = terrainTile->getMesh();
+  WrapperMesh mesh(mGridBounds, tileMesh, tileSizeX, tileSizeY);
+  heightfield.generateMesh(mesh, 0);
+  heightfield.clear();
+
+  // If we are not at the maximum zoom level we need to set child flags on the
+  // tile where child tiles overlap the dataset bounds.
+  if (coord.zoom != maxZoomLevel()) {
+    CRSBounds tileBounds = mGrid.tileBounds(coord);
+
+    if (! (bounds().overlaps(tileBounds))) {
+      terrainTile->setAllChildren(false);
+    } else {
+      if (bounds().overlaps(tileBounds.getSW())) {
+        terrainTile->setChildSW();
+      }
+      if (bounds().overlaps(tileBounds.getNW())) {
+        terrainTile->setChildNW();
+      }
+      if (bounds().overlaps(tileBounds.getNE())) {
+        terrainTile->setChildNE();
+      }
+      if (bounds().overlaps(tileBounds.getSE())) {
+        terrainTile->setChildSE();
+      }
+    }
+  }
+}
+
+MeshTile *
+ctb::MeshTiler::createMesh(GDALDataset *dataset, const TileCoordinate &coord) const {
+  // Copy the raster data into an array
+  float *rasterHeights = ctb::GDALDatasetReader::readRasterHeights(*this, dataset, coord, mGrid.tileSize(), mGrid.tileSize());
+
+  // Get a mesh tile represented by the tile coordinate
+  MeshTile *terrainTile = new MeshTile(coord);
+  prepareSettingsOfTile(terrainTile, dataset, coord, rasterHeights, mGrid.tileSize(), mGrid.tileSize());
+  CPLFree(rasterHeights);
+
+  return terrainTile;
+}
+
+MeshTile *
+ctb::MeshTiler::createMesh(GDALDataset *dataset, const TileCoordinate &coord, ctb::GDALDatasetReader *reader) const {
+  // Copy the raster data into an array
+  float *rasterHeights = reader->readRasterHeights(dataset, coord, mGrid.tileSize(), mGrid.tileSize());
+
+  // Get a mesh tile represented by the tile coordinate
+  MeshTile *terrainTile = new MeshTile(coord);
+  prepareSettingsOfTile(terrainTile, dataset, coord, rasterHeights, mGrid.tileSize(), mGrid.tileSize());
+  CPLFree(rasterHeights);
+
+  return terrainTile;
+}
+
+MeshTiler &
+ctb::MeshTiler::operator=(const MeshTiler &other) {
+  TerrainTiler::operator=(other);
+
+  return *this;
+}
+
+double ctb::MeshTiler::getEstimatedLevelZeroGeometricErrorForAHeightmap(
+  double maximumRadius,
+  double heightmapTerrainQuality,
+  int tileWidth,
+  int numberOfTilesAtLevelZero)
+{
+  double error = maximumRadius * 2 * M_PI * heightmapTerrainQuality;
+  error /= (double)(tileWidth * numberOfTilesAtLevelZero);
+  return error;
+}
diff --git a/src/MeshTiler.hpp b/src/MeshTiler.hpp
new file mode 100644
index 0000000..8207e25
--- /dev/null
+++ b/src/MeshTiler.hpp
@@ -0,0 +1,87 @@
+#ifndef MESHTILER_HPP
+#define MESHTILER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file MeshTiler.hpp
+ * @brief This declares the `MeshTiler` class
+ * @author Alvaro Huarte <ahuarte47@yahoo.es>
+ */
+
+#include "MeshTile.hpp"
+#include "TerrainTiler.hpp"
+
+namespace ctb {
+  class MeshTiler;
+}
+
+/**
+ * @brief Create `MeshTile`s from a GDAL Dataset
+ *
+ * This class derives from `GDALTiler` and `TerrainTiler` enabling `MeshTile`s 
+ * to be created for a specific `TileCoordinate`.
+ */
+class CTB_DLL ctb::MeshTiler :
+  public TerrainTiler
+{
+public:
+
+  /// Instantiate a tiler with all required arguments
+  MeshTiler(GDALDataset *poDataset, const Grid &grid, const TilerOptions &options, double meshQualityFactor = 1.0):
+    TerrainTiler(poDataset, grid, options),
+    mMeshQualityFactor(meshQualityFactor) {}
+
+  /// Instantiate a tiler with an empty GDAL dataset
+  MeshTiler(double meshQualityFactor = 1.0):
+    TerrainTiler(),
+    mMeshQualityFactor(meshQualityFactor) {}
+
+  /// Instantiate a tiler with a dataset and grid but no options
+  MeshTiler(GDALDataset *poDataset, const Grid &grid, double meshQualityFactor = 1.0):
+    TerrainTiler(poDataset, grid, TilerOptions()),
+    mMeshQualityFactor(meshQualityFactor) {}
+
+  /// Overload the assignment operator
+  MeshTiler &
+  operator=(const MeshTiler &other);
+
+  /// Create a mesh from a tile coordinate
+  MeshTile *
+  createMesh(GDALDataset *dataset, const TileCoordinate &coord) const;
+
+  /// Create a mesh from a tile coordinate
+  MeshTile *
+  createMesh(GDALDataset *dataset, const TileCoordinate &coord, GDALDatasetReader *reader) const;
+
+protected:
+
+  // Specifies the factor of the quality to convert terrain heightmaps to meshes.
+  double mMeshQualityFactor;
+
+  // Determines an appropriate geometric error estimate when the geometry comes from a heightmap.
+  static double getEstimatedLevelZeroGeometricErrorForAHeightmap(
+    double maximumRadius, 
+    double heightmapTerrainQuality, 
+    int tileWidth, 
+    int numberOfTilesAtLevelZero);
+
+  /// Assigns settings of Tile just to use.
+  void prepareSettingsOfTile(MeshTile *tile, GDALDataset *dataset, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const;
+};
+
+#endif /* MESHTILER_HPP */
diff --git a/src/RasterTiler.hpp b/src/RasterTiler.hpp
index f673e33..2bab229 100644
--- a/src/RasterTiler.hpp
+++ b/src/RasterTiler.hpp
@@ -54,8 +54,8 @@ class ctb::RasterTiler :
 
   /// Override to return a covariant data type
   virtual GDALTile *
-  createTile(const TileCoordinate &coord) const override {
-    return createRasterTile(coord);
+  createTile(GDALDataset *dataset, const TileCoordinate &coord) const override {
+    return createRasterTile(dataset, coord);
   }
 };
 
diff --git a/src/TerrainIterator.hpp b/src/TerrainIterator.hpp
index f009e00..8d8396b 100644
--- a/src/TerrainIterator.hpp
+++ b/src/TerrainIterator.hpp
@@ -56,6 +56,10 @@ class ctb::TerrainIterator :
   operator*() const override {
     return static_cast<TerrainTile *>(TilerIterator::operator*());
   }
+  virtual TerrainTile *
+  operator*(ctb::GDALDatasetReader *reader) const {
+    return (static_cast<const TerrainTiler &>(tiler)).createTile(tiler.dataset(), *(GridIterator::operator*()), reader);
+  }
 };
 
 #endif /* TERRAINITERATOR_HPP */
diff --git a/src/TerrainSerializer.hpp b/src/TerrainSerializer.hpp
new file mode 100644
index 0000000..521f6a3
--- /dev/null
+++ b/src/TerrainSerializer.hpp
@@ -0,0 +1,50 @@
+#ifndef TERRAINSERIALIZER_HPP
+#define TERRAINSERIALIZER_HPP
+
+/*******************************************************************************
+ * Copyright 2018 GeoData <geodata@soton.ac.uk>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may not
+ * use this file except in compliance with the License.  You may obtain a copy
+ * of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS, WITHOUT
+ * WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.  See the
+ * License for the specific language governing permissions and limitations
+ * under the License.
+ *******************************************************************************/
+
+/**
+ * @file TerrainSerializer.hpp
+ * @brief This declares and defines the `TerrainSerializer` class
+ */
+
+#include "config.hpp"
+#include "TileCoordinate.hpp"
+#include "TerrainTile.hpp"
+
+namespace ctb {
+  class TerrainSerializer;
+}
+
+/// Store `TerrainTile`s from a GDAL Dataset
+class CTB_DLL ctb::TerrainSerializer {
+public:
+
+  /// Start a new serialization task
+  virtual void startSerialization() = 0;
+
+  /// Returns if the specified Tile Coordinate should be serialized
+  virtual bool mustSerializeCoordinate(const ctb::TileCoordinate *coordinate) = 0;
+
+  /// Serialize a TerrainTile to the store
+  virtual bool serializeTile(const ctb::TerrainTile *tile) = 0;
+
+  /// Serialization finished, releases any resources loaded
+  virtual void endSerialization() = 0;
+};
+
+#endif /* TERRAINSERIALIZER_HPP */
diff --git a/src/TerrainTile.cpp b/src/TerrainTile.cpp
index 7b78262..7c416bf 100644
--- a/src/TerrainTile.cpp
+++ b/src/TerrainTile.cpp
@@ -28,6 +28,8 @@
 #include "TerrainTile.hpp"
 #include "GlobalGeodetic.hpp"
 #include "Bounds.hpp"
+#include "CTBFileOutputStream.hpp"
+#include "CTBZOutputStream.hpp"
 
 using namespace ctb;
 
@@ -136,9 +138,8 @@ Terrain::readFile(const char *fileName) {
  */
 void
 Terrain::writeFile(FILE *fp) const {
-  fwrite(mHeights.data(), TILE_CELL_SIZE * 2, 1, fp);
-  fwrite(&mChildren, 1, 1, fp);
-  fwrite(mMask, mMaskLength, 1, fp);
+  CTBFileOutputStream ostream(fp);
+  writeFile(ostream);
 }
 
 /**
@@ -146,41 +147,30 @@ Terrain::writeFile(FILE *fp) const {
  */
 void 
 Terrain::writeFile(const char *fileName) const {
-  gzFile terrainFile = gzopen(fileName, "wb");
+  CTBZFileOutputStream ostream(fileName);
+  writeFile(ostream);
+}
 
-  if (terrainFile == NULL) {
-    throw CTBException("Failed to open file");
-  }
+/**
+ * @details This writes raw terrain data to an output stream.
+ */
+void
+Terrain::writeFile(CTBOutputStream &ostream) const {
 
   // Write the height data
-  if (gzwrite(terrainFile, mHeights.data(), TILE_CELL_SIZE * 2) == 0) {
-    gzclose(terrainFile);
+  if (ostream.write((const void *)mHeights.data(), TILE_CELL_SIZE * 2) != TILE_CELL_SIZE * 2) {
     throw CTBException("Failed to write height data");
   }
 
   // Write the child flags
-  if (gzputc(terrainFile, mChildren) == -1) {
-    gzclose(terrainFile);
+  if (ostream.write(&mChildren, 1) != 1) {
     throw CTBException("Failed to write child flags");
   }
 
   // Write the water mask
-  if (gzwrite(terrainFile, mMask, mMaskLength) == 0) {
-    gzclose(terrainFile);
+  if (ostream.write(&mMask, mMaskLength) != mMaskLength) {
     throw CTBException("Failed to write water mask");
   }
-
-  // Try and close the file
-  switch (gzclose(terrainFile)) {
-  case Z_OK:
-    break;
-  case Z_STREAM_ERROR:
-  case Z_ERRNO:
-  case Z_MEM_ERROR:
-  case Z_BUF_ERROR:
-  default:
-    throw CTBException("Failed to close file");
-  }
 }
 
 std::vector<bool>
@@ -334,6 +324,11 @@ TerrainTile::heightsToRaster() const {
 
   // Create the spatial reference system for the raster
   OGRSpatialReference oSRS;
+
+  #if ( GDAL_VERSION_MAJOR >= 3 )
+  oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+  #endif
+
   if (oSRS.importFromEPSG(4326) != OGRERR_NONE) {
     throw CTBException("Could not create EPSG:4326 spatial reference");
   }
diff --git a/src/TerrainTile.hpp b/src/TerrainTile.hpp
index 5941166..ce739ad 100644
--- a/src/TerrainTile.hpp
+++ b/src/TerrainTile.hpp
@@ -29,6 +29,7 @@
 #include "config.hpp"
 #include "Tile.hpp"
 #include "TileCoordinate.hpp"
+#include "CTBOutputStream.hpp"
 
 namespace ctb {
   class Terrain;
@@ -65,6 +66,10 @@ class CTB_DLL ctb::Terrain {
   void
   writeFile(const char *fileName) const;
 
+  /// Write terrain data to an output stream
+  void
+  writeFile(CTBOutputStream &ostream) const;
+
   /// Get the water mask as a boolean mask
   std::vector<bool>
   mask();
diff --git a/src/TerrainTiler.cpp b/src/TerrainTiler.cpp
index 174abf0..77463aa 100644
--- a/src/TerrainTiler.cpp
+++ b/src/TerrainTiler.cpp
@@ -21,32 +21,20 @@
 
 #include "CTBException.hpp"
 #include "TerrainTiler.hpp"
+#include "GDALDatasetReader.hpp"
 
 using namespace ctb;
 
-TerrainTile *
-ctb::TerrainTiler::createTile(const TileCoordinate &coord) const {
-  // Get a terrain tile represented by the tile coordinate
-  TerrainTile *terrainTile = new TerrainTile(coord);
-  GDALTile *rasterTile = createRasterTile(coord); // the raster associated with this tile coordinate
-  GDALRasterBand *heightsBand = rasterTile->dataset->GetRasterBand(1);
-
-  // Copy the raster data into an array
-  float rasterHeights[TerrainTile::TILE_CELL_SIZE];
-  if (heightsBand->RasterIO(GF_Read, 0, 0, TILE_SIZE, TILE_SIZE,
-                            (void *) rasterHeights, TILE_SIZE, TILE_SIZE, GDT_Float32,
-                            0, 0) != CE_None) {
-    throw CTBException("Could not read heights from raster");
-  }
-
-  delete rasterTile;
+void 
+ctb::TerrainTiler::prepareSettingsOfTile(TerrainTile *terrainTile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const {
+  const ctb::i_tile TILE_CELL_SIZE = tileSizeX * tileSizeY;
 
   // Convert the raster data into the terrain tile heights.  This assumes the
   // input raster data represents meters above sea level. Each terrain height
   // value is the number of 1/5 meter units above -1000 meters.
   // TODO: try doing this using a VRT derived band:
   // (http://www.gdal.org/gdal_vrttut.html)
-  for (unsigned short int i = 0; i < TerrainTile::TILE_CELL_SIZE; i++) {
+  for (unsigned short int i = 0; i < TILE_CELL_SIZE; i++) {
     terrainTile->mHeights[i] = (i_terrain_height) ((rasterHeights[i] + 1000) * 5);
   }
 
@@ -72,14 +60,38 @@ ctb::TerrainTiler::createTile(const TileCoordinate &coord) const {
       }
     }
   }
+}
+
+TerrainTile *
+ctb::TerrainTiler::createTile(GDALDataset *dataset, const TileCoordinate &coord) const {
+  // Copy the raster data into an array
+  float *rasterHeights = ctb::GDALDatasetReader::readRasterHeights(*this, dataset, coord, TILE_SIZE, TILE_SIZE);
+
+  // Get a terrain tile represented by the tile coordinate
+  TerrainTile *terrainTile = new TerrainTile(coord);
+  prepareSettingsOfTile(terrainTile, coord, rasterHeights, TILE_SIZE, TILE_SIZE);
+  CPLFree(rasterHeights);
+
+  return terrainTile;
+}
+
+TerrainTile *
+ctb::TerrainTiler::createTile(GDALDataset *dataset, const TileCoordinate &coord, ctb::GDALDatasetReader *reader) const {
+  // Copy the raster data into an array
+  float *rasterHeights = reader->readRasterHeights(dataset, coord, TILE_SIZE, TILE_SIZE);
+
+  // Get a mesh tile represented by the tile coordinate
+  TerrainTile *terrainTile = new TerrainTile(coord);
+  prepareSettingsOfTile(terrainTile, coord, rasterHeights, TILE_SIZE, TILE_SIZE);
+  CPLFree(rasterHeights);
 
   return terrainTile;
 }
 
 GDALTile *
-ctb::TerrainTiler::createRasterTile(const TileCoordinate &coord) const {
+ctb::TerrainTiler::createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const {
   // Ensure we have some data from which to create a tile
-  if (poDataset && poDataset->GetRasterCount() < 1) {
+  if (dataset && dataset->GetRasterCount() < 1) {
     throw CTBException("At least one band must be present in the GDAL dataset");
   }
 
@@ -97,7 +109,7 @@ ctb::TerrainTiler::createRasterTile(const TileCoordinate &coord) const {
   adfGeoTransform[4] = 0;
   adfGeoTransform[5] = -resolution;
 
-  GDALTile *tile = GDALTiler::createRasterTile(adfGeoTransform);
+  GDALTile *tile = GDALTiler::createRasterTile(dataset, adfGeoTransform);
 
   // The previous geotransform represented the data with an overlap as required
   // by the terrain specification.  This now needs to be overwritten so that
diff --git a/src/TerrainTiler.hpp b/src/TerrainTiler.hpp
index 00da46b..171ae3e 100644
--- a/src/TerrainTiler.hpp
+++ b/src/TerrainTiler.hpp
@@ -59,13 +59,17 @@ class CTB_DLL ctb::TerrainTiler :
 
   /// Override to return a covariant data type
   TerrainTile *
-  createTile(const TileCoordinate &coord) const override;
+  createTile(GDALDataset *dataset, const TileCoordinate &coord) const override;
+
+  /// Create a tile from a tile coordinate
+  TerrainTile *
+  createTile(GDALDataset *dataset, const TileCoordinate &coord, GDALDatasetReader *reader) const;
 
 protected:
 
   /// Create a `GDALTile` representing the required terrain tile data
   virtual GDALTile *
-  createRasterTile(const TileCoordinate &coord) const override;
+  createRasterTile(GDALDataset *dataset, const TileCoordinate &coord) const override;
 
   /**
    * @brief Get terrain bounds shifted to introduce a pixel overlap
@@ -97,6 +101,9 @@ class CTB_DLL ctb::TerrainTiler :
 
     return tile;
   }
+
+  /// Assigns settings of Tile just to use.
+  void prepareSettingsOfTile(TerrainTile *tile, const TileCoordinate &coord, float *rasterHeights, ctb::i_tile tileSizeX, ctb::i_tile tileSizeY) const;
 };
 
 #endif /* TERRAINTILER_HPP */
diff --git a/src/TilerIterator.hpp b/src/TilerIterator.hpp
index fca6f3a..83ed2d9 100644
--- a/src/TilerIterator.hpp
+++ b/src/TilerIterator.hpp
@@ -55,7 +55,7 @@ class ctb::TilerIterator :
   /// Override the dereference operator to return a Tile
   virtual Tile *
   operator*() const {
-    return tiler.createTile(*(GridIterator::operator*()));
+    return tiler.createTile(tiler.dataset(), *(GridIterator::operator*()));
   }
 
 protected:
diff --git a/src/gdaloverviewdataset-gdal2x.cpp b/src/gdaloverviewdataset-gdal2x.cpp
new file mode 100644
index 0000000..bfb4ca5
--- /dev/null
+++ b/src/gdaloverviewdataset-gdal2x.cpp
@@ -0,0 +1,618 @@
+/******************************************************************************
+ *
+ * Project:  GDAL Core
+ * Purpose:  Implementation of a dataset overview warping class
+ * Author:   Even Rouault, <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2014, Even Rouault, <even dot rouault at spatialys dot com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include "cpl_port.h"
+#include "gdal_priv.h"
+
+#include <cstring>
+
+#include "cpl_conv.h"
+#include "cpl_error.h"
+#include "cpl_progress.h"
+#include "cpl_string.h"
+#include "gdal.h"
+#include "gdal_mdreader.h"
+#include "gdal_proxy.h"
+
+CPL_CVSID("$Id$")
+
+/** In GDAL, GDALRasterBand::GetOverview() returns a stand-alone band, that may
+    have no parent dataset. This can be inconvenient in certain contexts, where
+    cross-band processing must be done, or when API expect a fully fledged
+    dataset.  Furthermore even if overview band has a container dataset, that
+    one often fails to declare its projection, geotransform, etc... which make
+    it somehow useless. GDALOverviewDataset remedies to those deficiencies.
+*/
+
+class GDALOverviewBand;
+
+/* ******************************************************************** */
+/*                          GDALOverviewDataset                         */
+/* ******************************************************************** */
+
+class GDALOverviewDataset : public GDALDataset
+{
+  private:
+    friend class GDALOverviewBand;
+
+    GDALDataset* poMainDS;
+
+    GDALDataset* poOvrDS;  // Will be often NULL.
+    int          nOvrLevel;
+    int          bThisLevelOnly;
+
+    int          nGCPCount;
+    GDAL_GCP    *pasGCPList;
+    char       **papszMD_RPC;
+    char       **papszMD_GEOLOCATION;
+
+    static void  Rescale( char**& papszMD, const char* pszItem,
+                          double dfRatio, double dfDefaultVal );
+
+  protected:
+    CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
+                      void *, int, int, GDALDataType,
+                      int, int *,
+                      GSpacing, GSpacing, GSpacing,
+                      GDALRasterIOExtraArg* psExtraArg ) override;
+
+  public:
+    GDALOverviewDataset( GDALDataset* poMainDS,
+                         int nOvrLevel,
+                         int bThisLevelOnly );
+    ~GDALOverviewDataset() override;
+
+    const char *GetProjectionRef( void ) override;
+    CPLErr GetGeoTransform( double * ) override;
+
+    int GetGCPCount() override;
+    const char *GetGCPProjection() override;
+    const GDAL_GCP *GetGCPs() override;
+
+    char  **GetMetadata( const char * pszDomain = "" ) override;
+    const char *GetMetadataItem( const char * pszName,
+                                 const char * pszDomain = "" ) override;
+
+    int CloseDependentDatasets() override;
+
+  private:
+    CPL_DISALLOW_COPY_ASSIGN(GDALOverviewDataset)
+};
+
+/* ******************************************************************** */
+/*                           GDALOverviewBand                           */
+/* ******************************************************************** */
+
+class GDALOverviewBand : public GDALProxyRasterBand
+{
+  protected:
+    friend class GDALOverviewDataset;
+
+    GDALRasterBand*         poUnderlyingBand;
+    GDALRasterBand* RefUnderlyingRasterBand() override;
+
+  public:
+    GDALOverviewBand( GDALOverviewDataset* poDS, int nBand );
+    ~GDALOverviewBand() override;
+
+    CPLErr FlushCache() override;
+
+    int GetOverviewCount() override;
+    GDALRasterBand *GetOverview( int ) override;
+
+  private:
+    CPL_DISALLOW_COPY_ASSIGN(GDALOverviewBand)
+};
+
+/************************************************************************/
+/*                       GDALCreateOverviewDataset()                    */
+/************************************************************************/
+
+// Takes a reference on poMainDS in case of success.
+GDALDataset* GDALCreateOverviewDataset( GDALDataset* poMainDS, int nOvrLevel,
+                                        int bThisLevelOnly )
+{
+    // Sanity checks.
+    const int nBands = poMainDS->GetRasterCount();
+    if( nBands == 0 )
+        return nullptr;
+
+    for( int i = 1; i<= nBands; ++i )
+    {
+        if( poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel) == nullptr )
+        {
+            return nullptr;
+        }
+        if( poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel)->GetXSize() !=
+            poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetXSize() ||
+            poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel)->GetYSize() !=
+            poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetYSize() )
+        {
+            return nullptr;
+        }
+    }
+
+    return new GDALOverviewDataset(poMainDS, nOvrLevel, bThisLevelOnly);
+}
+
+/************************************************************************/
+/*                        GDALOverviewDataset()                         */
+/************************************************************************/
+
+GDALOverviewDataset::GDALOverviewDataset( GDALDataset* poMainDSIn,
+                                          int nOvrLevelIn,
+                                          int bThisLevelOnlyIn ) :
+    poMainDS(poMainDSIn),
+    nOvrLevel(nOvrLevelIn),
+    bThisLevelOnly(bThisLevelOnlyIn),
+    nGCPCount(0),
+    pasGCPList(nullptr),
+    papszMD_RPC(nullptr),
+    papszMD_GEOLOCATION(nullptr)
+{
+    poMainDSIn->Reference();
+    eAccess = poMainDS->GetAccess();
+    nRasterXSize =
+        poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetXSize();
+    nRasterYSize =
+        poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetYSize();
+    poOvrDS = poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetDataset();
+    if( poOvrDS != nullptr && poOvrDS == poMainDS )
+    {
+        CPLDebug( "GDAL",
+                  "Dataset of overview is the same as the main band. "
+                  "This is not expected");
+        poOvrDS = nullptr;
+    }
+    nBands = poMainDS->GetRasterCount();
+    for( int i = 0; i < nBands; ++i )
+    {
+        SetBand(i+1, new GDALOverviewBand(this, i+1));
+    }
+
+    // We create a fake driver that has the same name as the original
+    // one, but we cannot use the real driver object, so that code
+    // doesn't try to cast the GDALOverviewDataset* as a native dataset
+    // object.
+    if( poMainDS->GetDriver() != nullptr )
+    {
+        poDriver = new GDALDriver();
+        poDriver->SetDescription(poMainDS->GetDriver()->GetDescription());
+        poDriver->SetMetadata(poMainDS->GetDriver()->GetMetadata());
+    }
+
+    SetDescription( poMainDS->GetDescription() );
+
+    CPLDebug( "GDAL", "GDALOverviewDataset(%s, this=%p) creation.",
+              poMainDS->GetDescription(), this );
+
+    papszOpenOptions = CSLDuplicate(poMainDS->GetOpenOptions());
+    // Add OVERVIEW_LEVEL if not called from GDALOpenEx(), but directly.
+    papszOpenOptions = CSLSetNameValue(papszOpenOptions, "OVERVIEW_LEVEL",
+                                       CPLSPrintf("%d", nOvrLevel));
+}
+
+/************************************************************************/
+/*                       ~GDALOverviewDataset()                         */
+/************************************************************************/
+
+GDALOverviewDataset::~GDALOverviewDataset()
+{
+    FlushCache();
+
+    CloseDependentDatasets();
+
+    if( nGCPCount > 0 )
+    {
+        GDALDeinitGCPs( nGCPCount, pasGCPList );
+        CPLFree( pasGCPList );
+    }
+    CSLDestroy(papszMD_RPC);
+
+    CSLDestroy(papszMD_GEOLOCATION);
+
+    delete poDriver;
+}
+
+/************************************************************************/
+/*                      CloseDependentDatasets()                        */
+/************************************************************************/
+
+int GDALOverviewDataset::CloseDependentDatasets()
+{
+    bool bRet = false;
+
+    if( poMainDS )
+    {
+        for( int i = 0; i < nBands; ++i )
+        {
+            GDALOverviewBand* const band =
+                dynamic_cast<GDALOverviewBand*>(papoBands[i]);
+            if( band == nullptr )
+            {
+                CPLError( CE_Fatal, CPLE_AppDefined,
+                            "OverviewBand cast fail." );
+                return false;
+            }
+            band->poUnderlyingBand = nullptr;
+        }
+        if( poMainDS->ReleaseRef() )
+            bRet = true;
+        poMainDS = nullptr;
+    }
+
+    return bRet;
+}
+
+/************************************************************************/
+/*                             IRasterIO()                              */
+/*                                                                      */
+/*      The default implementation of IRasterIO() is to pass the        */
+/*      request off to each band objects rasterio methods with          */
+/*      appropriate arguments.                                          */
+/************************************************************************/
+
+CPLErr GDALOverviewDataset::IRasterIO( GDALRWFlag eRWFlag,
+                                       int nXOff, int nYOff,
+                                       int nXSize, int nYSize,
+                                       void * pData,
+                                       int nBufXSize, int nBufYSize,
+                                       GDALDataType eBufType,
+                                       int nBandCount, int *panBandMap,
+                                       GSpacing nPixelSpace,
+                                       GSpacing nLineSpace,
+                                       GSpacing nBandSpace,
+                                       GDALRasterIOExtraArg* psExtraArg )
+
+{
+    // In case the overview bands are really linked to a dataset, then issue
+    // the request to that dataset.
+    if( poOvrDS != nullptr )
+    {
+        return poOvrDS->RasterIO(
+            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
+            eBufType, nBandCount, panBandMap, nPixelSpace,
+            nLineSpace, nBandSpace,
+            psExtraArg);
+    }
+
+    GDALProgressFunc  pfnProgressGlobal = psExtraArg->pfnProgress;
+    void *pProgressDataGlobal = psExtraArg->pProgressData;
+    CPLErr eErr = CE_None;
+
+    for( int iBandIndex = 0;
+         iBandIndex < nBandCount && eErr == CE_None;
+         ++iBandIndex )
+    {
+        GDALOverviewBand *poBand =
+            dynamic_cast<GDALOverviewBand *>(
+                GetRasterBand(panBandMap[iBandIndex]) );
+        if( poBand == nullptr )
+        {
+            eErr = CE_Failure;
+            break;
+        }
+
+        GByte *pabyBandData =
+            static_cast<GByte *>(pData) + iBandIndex * nBandSpace;
+
+        psExtraArg->pfnProgress = GDALScaledProgress;
+        psExtraArg->pProgressData =
+            GDALCreateScaledProgress( 1.0 * iBandIndex / nBandCount,
+                                      1.0 * (iBandIndex + 1) / nBandCount,
+                                      pfnProgressGlobal,
+                                      pProgressDataGlobal );
+
+        eErr = poBand->IRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize,
+                                  pabyBandData,
+                                  nBufXSize, nBufYSize,
+                                  eBufType, nPixelSpace,
+                                  nLineSpace, psExtraArg );
+
+        GDALDestroyScaledProgress( psExtraArg->pProgressData );
+    }
+
+    psExtraArg->pfnProgress = pfnProgressGlobal;
+    psExtraArg->pProgressData = pProgressDataGlobal;
+
+    return eErr;
+}
+
+/************************************************************************/
+/*                          GetProjectionRef()                          */
+/************************************************************************/
+
+const char *GDALOverviewDataset::GetProjectionRef()
+
+{
+    return poMainDS->GetProjectionRef();
+}
+
+/************************************************************************/
+/*                          GetGeoTransform()                           */
+/************************************************************************/
+
+CPLErr GDALOverviewDataset::GetGeoTransform( double * padfTransform )
+
+{
+    double adfGeoTransform[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+    if( poMainDS->GetGeoTransform(adfGeoTransform) != CE_None )
+        return CE_Failure;
+
+    adfGeoTransform[1] *=
+        static_cast<double>(poMainDS->GetRasterXSize()) / nRasterXSize;
+    adfGeoTransform[2] *=
+        static_cast<double>(poMainDS->GetRasterYSize()) / nRasterYSize;
+    adfGeoTransform[4] *=
+        static_cast<double>(poMainDS->GetRasterXSize()) / nRasterXSize;
+    adfGeoTransform[5] *=
+        static_cast<double>(poMainDS->GetRasterYSize()) / nRasterYSize;
+
+    memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
+
+    return CE_None;
+}
+
+/************************************************************************/
+/*                            GetGCPCount()                             */
+/************************************************************************/
+
+int GDALOverviewDataset::GetGCPCount()
+
+{
+    return poMainDS->GetGCPCount();
+}
+
+/************************************************************************/
+/*                          GetGCPProjection()                          */
+/************************************************************************/
+
+const char *GDALOverviewDataset::GetGCPProjection()
+
+{
+    return poMainDS->GetGCPProjection();
+}
+
+/************************************************************************/
+/*                               GetGCPs()                              */
+/************************************************************************/
+
+const GDAL_GCP *GDALOverviewDataset::GetGCPs()
+
+{
+    if( pasGCPList != nullptr )
+        return pasGCPList;
+
+    const GDAL_GCP* pasGCPsMain = poMainDS->GetGCPs();
+    if( pasGCPsMain == nullptr )
+        return nullptr;
+    nGCPCount = poMainDS->GetGCPCount();
+
+    pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPsMain );
+    for( int i = 0; i < nGCPCount; ++i )
+    {
+        pasGCPList[i].dfGCPPixel *= static_cast<double>(nRasterXSize) /
+            poMainDS->GetRasterXSize();
+        pasGCPList[i].dfGCPLine *= static_cast<double>(nRasterYSize) /
+            poMainDS->GetRasterYSize();
+    }
+    return pasGCPList;
+}
+
+/************************************************************************/
+/*                             Rescale()                                */
+/************************************************************************/
+
+void GDALOverviewDataset::Rescale( char**& papszMD, const char* pszItem,
+                                   double dfRatio, double dfDefaultVal )
+{
+    double dfVal =
+        CPLAtofM( CSLFetchNameValueDef(papszMD, pszItem,
+                                       CPLSPrintf("%.18g", dfDefaultVal)) );
+    dfVal *= dfRatio;
+    papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.18g", dfVal));
+}
+
+/************************************************************************/
+/*                            GetMetadata()                             */
+/************************************************************************/
+
+char  **GDALOverviewDataset::GetMetadata( const char * pszDomain )
+{
+    if( poOvrDS != nullptr )
+    {
+        char** papszMD = poOvrDS->GetMetadata(pszDomain);
+        if( papszMD != nullptr )
+            return papszMD;
+    }
+
+    char** papszMD = poMainDS->GetMetadata(pszDomain);
+
+    // We may need to rescale some values from the RPC metadata domain.
+    if( pszDomain != nullptr && EQUAL(pszDomain, MD_DOMAIN_RPC) &&
+        papszMD != nullptr )
+    {
+        if( papszMD_RPC )
+            return papszMD_RPC;
+        papszMD_RPC = CSLDuplicate(papszMD);
+
+        Rescale( papszMD_RPC, RPC_LINE_OFF,
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 0.0 );
+        Rescale( papszMD_RPC, RPC_LINE_SCALE,
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 1.0 );
+        Rescale( papszMD_RPC, RPC_SAMP_OFF,
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 0.0 );
+        Rescale( papszMD_RPC, RPC_SAMP_SCALE,
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 1.0 );
+
+        papszMD = papszMD_RPC;
+    }
+
+    // We may need to rescale some values from the GEOLOCATION metadata domain.
+    if( pszDomain != nullptr && EQUAL(pszDomain, "GEOLOCATION") &&
+        papszMD != nullptr )
+    {
+        if( papszMD_GEOLOCATION )
+            return papszMD_GEOLOCATION;
+        papszMD_GEOLOCATION = CSLDuplicate(papszMD);
+
+        Rescale( papszMD_GEOLOCATION, "PIXEL_OFFSET",
+                 static_cast<double>(poMainDS->GetRasterXSize()) /
+                 nRasterXSize, 0.0 );
+        Rescale( papszMD_GEOLOCATION, "LINE_OFFSET",
+                 static_cast<double>(poMainDS->GetRasterYSize()) /
+                 nRasterYSize, 0.0 );
+
+        Rescale( papszMD_GEOLOCATION, "PIXEL_STEP",
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 1.0 );
+        Rescale( papszMD_GEOLOCATION, "LINE_STEP",
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 1.0 );
+
+        papszMD = papszMD_GEOLOCATION;
+    }
+
+    return papszMD;
+}
+
+/************************************************************************/
+/*                          GetMetadataItem()                           */
+/************************************************************************/
+
+const char *GDALOverviewDataset::GetMetadataItem( const char * pszName,
+                                                  const char * pszDomain )
+{
+    if( poOvrDS != nullptr )
+    {
+        const char* pszValue = poOvrDS->GetMetadataItem(pszName, pszDomain);
+        if( pszValue != nullptr )
+            return pszValue;
+    }
+
+    if( pszDomain != nullptr && (EQUAL(pszDomain, "RPC") ||
+                              EQUAL(pszDomain, "GEOLOCATION")) )
+    {
+        char** papszMD = GetMetadata(pszDomain);
+        return CSLFetchNameValue(papszMD, pszName);
+    }
+
+    return poMainDS->GetMetadataItem(pszName, pszDomain);
+}
+
+/************************************************************************/
+/*                          GDALOverviewBand()                          */
+/************************************************************************/
+
+GDALOverviewBand::GDALOverviewBand( GDALOverviewDataset* poDSIn, int nBandIn ) :
+    poUnderlyingBand(poDSIn->poMainDS->GetRasterBand(nBandIn)->
+                         GetOverview(poDSIn->nOvrLevel))
+{
+    poDS = poDSIn;
+    nBand = nBandIn;
+    nRasterXSize = poDSIn->nRasterXSize;
+    nRasterYSize = poDSIn->nRasterYSize;
+    eDataType = poUnderlyingBand->GetRasterDataType();
+    poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
+}
+
+/************************************************************************/
+/*                         ~GDALOverviewBand()                          */
+/************************************************************************/
+
+GDALOverviewBand::~GDALOverviewBand()
+{
+    FlushCache();
+}
+
+/************************************************************************/
+/*                              FlushCache()                            */
+/************************************************************************/
+
+CPLErr GDALOverviewBand::FlushCache()
+{
+    if( poUnderlyingBand )
+        return poUnderlyingBand->FlushCache();
+    return CE_None;
+}
+
+/************************************************************************/
+/*                        RefUnderlyingRasterBand()                     */
+/************************************************************************/
+
+GDALRasterBand* GDALOverviewBand::RefUnderlyingRasterBand()
+{
+    if( poUnderlyingBand )
+        return poUnderlyingBand;
+
+    return nullptr;
+}
+
+/************************************************************************/
+/*                         GetOverviewCount()                           */
+/************************************************************************/
+
+int GDALOverviewBand::GetOverviewCount()
+{
+    GDALOverviewDataset * const poOvrDS =
+        dynamic_cast<GDALOverviewDataset *>(poDS);
+    if( poOvrDS == nullptr )
+    {
+        CPLError( CE_Fatal, CPLE_AppDefined, "OverviewDataset cast fail." );
+        return 0;
+    }
+    if( poOvrDS->bThisLevelOnly )
+        return 0;
+    GDALDataset * const poMainDS = poOvrDS->poMainDS;
+    return poMainDS->GetRasterBand(nBand)->GetOverviewCount()
+        - poOvrDS->nOvrLevel - 1;
+}
+
+/************************************************************************/
+/*                           GetOverview()                              */
+/************************************************************************/
+
+GDALRasterBand *GDALOverviewBand::GetOverview( int iOvr )
+{
+    if( iOvr < 0 || iOvr >= GetOverviewCount() )
+        return nullptr;
+    GDALOverviewDataset * const poOvrDS =
+        dynamic_cast<GDALOverviewDataset *>(poDS);
+    if( poOvrDS == nullptr )
+    {
+        CPLError( CE_Fatal, CPLE_AppDefined, "OverviewDataset cast fail." );
+        return nullptr;
+    }
+    GDALDataset * const poMainDS = poOvrDS->poMainDS;
+    return poMainDS->GetRasterBand(nBand)->
+        GetOverview(iOvr + poOvrDS->nOvrLevel + 1);
+}
diff --git a/src/gdaloverviewdataset.cpp b/src/gdaloverviewdataset.cpp
new file mode 100644
index 0000000..783f59a
--- /dev/null
+++ b/src/gdaloverviewdataset.cpp
@@ -0,0 +1,671 @@
+/******************************************************************************
+ *
+ * Project:  GDAL Core
+ * Purpose:  Implementation of a dataset overview warping class
+ * Author:   Even Rouault, <even dot rouault at spatialys dot com>
+ *
+ ******************************************************************************
+ * Copyright (c) 2014, Even Rouault, <even dot rouault at spatialys dot com>
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ ****************************************************************************/
+
+#include "cpl_port.h"
+#include "gdal_priv.h"
+
+#include <cstring>
+
+#include "cpl_conv.h"
+#include "cpl_error.h"
+#include "cpl_progress.h"
+#include "cpl_string.h"
+#include "gdal.h"
+#include "gdal_mdreader.h"
+#include "gdal_proxy.h"
+
+CPL_CVSID("$Id$")
+
+/** In GDAL, GDALRasterBand::GetOverview() returns a stand-alone band, that may
+    have no parent dataset. This can be inconvenient in certain contexts, where
+    cross-band processing must be done, or when API expect a fully fledged
+    dataset.  Furthermore even if overview band has a container dataset, that
+    one often fails to declare its projection, geotransform, etc... which make
+    it somehow useless. GDALOverviewDataset remedies to those deficiencies.
+*/
+
+class GDALOverviewBand;
+
+/* ******************************************************************** */
+/*                          GDALOverviewDataset                         */
+/* ******************************************************************** */
+
+class GDALOverviewDataset final: public GDALDataset
+{
+  private:
+    friend class GDALOverviewBand;
+
+    GDALDataset* poMainDS = nullptr;
+
+    GDALDataset* poOvrDS = nullptr;  // Will be often NULL.
+    int          nOvrLevel = 0;
+    int          bThisLevelOnly = 0;
+
+    int          nGCPCount = 0;
+    GDAL_GCP    *pasGCPList = nullptr;
+    char       **papszMD_RPC = nullptr;
+    char       **papszMD_GEOLOCATION = nullptr;
+    GDALOverviewBand* m_poMaskBand = nullptr;
+
+    static void  Rescale( char**& papszMD, const char* pszItem,
+                          double dfRatio, double dfDefaultVal );
+
+  protected:
+    CPLErr IRasterIO( GDALRWFlag, int, int, int, int,
+                      void *, int, int, GDALDataType,
+                      int, int *,
+                      GSpacing, GSpacing, GSpacing,
+                      GDALRasterIOExtraArg* psExtraArg ) override;
+
+  public:
+    GDALOverviewDataset( GDALDataset* poMainDS,
+                         int nOvrLevel,
+                         int bThisLevelOnly );
+    ~GDALOverviewDataset() override;
+
+    const OGRSpatialReference* GetSpatialRef() const override;
+    CPLErr GetGeoTransform( double * ) override;
+
+    int GetGCPCount() override;
+    const OGRSpatialReference *GetGCPSpatialRef() const override;
+    const GDAL_GCP *GetGCPs() override;
+
+    char  **GetMetadata( const char * pszDomain = "" ) override;
+    const char *GetMetadataItem( const char * pszName,
+                                 const char * pszDomain = "" ) override;
+
+    int CloseDependentDatasets() override;
+
+  private:
+    CPL_DISALLOW_COPY_ASSIGN(GDALOverviewDataset)
+};
+
+/* ******************************************************************** */
+/*                           GDALOverviewBand                           */
+/* ******************************************************************** */
+
+class GDALOverviewBand final: public GDALProxyRasterBand
+{
+  protected:
+    friend class GDALOverviewDataset;
+
+    GDALRasterBand*         poUnderlyingBand = nullptr;
+    GDALRasterBand* RefUnderlyingRasterBand() override;
+
+  public:
+    GDALOverviewBand( GDALOverviewDataset* poDS, int nBand );
+    ~GDALOverviewBand() override;
+
+    CPLErr FlushCache() override;
+
+    int GetOverviewCount() override;
+    GDALRasterBand *GetOverview( int ) override;
+
+    int GetMaskFlags() override;
+    GDALRasterBand* GetMaskBand() override;
+
+  private:
+    CPL_DISALLOW_COPY_ASSIGN(GDALOverviewBand)
+};
+
+/************************************************************************/
+/*                       GDALCreateOverviewDataset()                    */
+/************************************************************************/
+
+// Takes a reference on poMainDS in case of success.
+GDALDataset* GDALCreateOverviewDataset( GDALDataset* poMainDS, int nOvrLevel,
+                                        int bThisLevelOnly )
+{
+    // Sanity checks.
+    const int nBands = poMainDS->GetRasterCount();
+    if( nBands == 0 )
+        return nullptr;
+
+    for( int i = 1; i<= nBands; ++i )
+    {
+        if( poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel) == nullptr )
+        {
+            return nullptr;
+        }
+        if( poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel)->GetXSize() !=
+            poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetXSize() ||
+            poMainDS->GetRasterBand(i)->GetOverview(nOvrLevel)->GetYSize() !=
+            poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetYSize() )
+        {
+            return nullptr;
+        }
+    }
+
+    return new GDALOverviewDataset(poMainDS, nOvrLevel, bThisLevelOnly);
+}
+
+/************************************************************************/
+/*                        GDALOverviewDataset()                         */
+/************************************************************************/
+
+GDALOverviewDataset::GDALOverviewDataset( GDALDataset* poMainDSIn,
+                                          int nOvrLevelIn,
+                                          int bThisLevelOnlyIn ) :
+    poMainDS(poMainDSIn),
+    nOvrLevel(nOvrLevelIn),
+    bThisLevelOnly(bThisLevelOnlyIn)
+{
+    poMainDSIn->Reference();
+    eAccess = poMainDS->GetAccess();
+    nRasterXSize =
+        poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetXSize();
+    nRasterYSize =
+        poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetYSize();
+    poOvrDS = poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetDataset();
+    if( poOvrDS != nullptr && poOvrDS == poMainDS )
+    {
+        CPLDebug( "GDAL",
+                  "Dataset of overview is the same as the main band. "
+                  "This is not expected");
+        poOvrDS = nullptr;
+    }
+    nBands = poMainDS->GetRasterCount();
+    for( int i = 0; i < nBands; ++i )
+    {
+        SetBand(i+1, new GDALOverviewBand(this, i+1));
+    }
+
+    if( poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetMaskFlags() == GMF_PER_DATASET )
+    {
+        auto poOvrMaskBand = poMainDS->GetRasterBand(1)->GetOverview(nOvrLevel)->GetMaskBand();
+        if( poOvrMaskBand  && poOvrMaskBand->GetXSize() == nRasterXSize &&
+            poOvrMaskBand->GetYSize() == nRasterYSize )
+        {
+            m_poMaskBand = new GDALOverviewBand(this, 0);
+        }
+    }
+
+    // We create a fake driver that has the same name as the original
+    // one, but we cannot use the real driver object, so that code
+    // doesn't try to cast the GDALOverviewDataset* as a native dataset
+    // object.
+    if( poMainDS->GetDriver() != nullptr )
+    {
+        poDriver = new GDALDriver();
+        poDriver->SetDescription(poMainDS->GetDriver()->GetDescription());
+        poDriver->SetMetadata(poMainDS->GetDriver()->GetMetadata());
+    }
+
+    SetDescription( poMainDS->GetDescription() );
+
+    CPLDebug( "GDAL", "GDALOverviewDataset(%s, this=%p) creation.",
+              poMainDS->GetDescription(), this );
+
+    papszOpenOptions = CSLDuplicate(poMainDS->GetOpenOptions());
+    // Add OVERVIEW_LEVEL if not called from GDALOpenEx(), but directly.
+    papszOpenOptions = CSLSetNameValue(papszOpenOptions, "OVERVIEW_LEVEL",
+                                       CPLSPrintf("%d", nOvrLevel));
+}
+
+/************************************************************************/
+/*                       ~GDALOverviewDataset()                         */
+/************************************************************************/
+
+GDALOverviewDataset::~GDALOverviewDataset()
+{
+    GDALOverviewDataset::FlushCache();
+
+    GDALOverviewDataset::CloseDependentDatasets();
+
+    if( nGCPCount > 0 )
+    {
+        GDALDeinitGCPs( nGCPCount, pasGCPList );
+        CPLFree( pasGCPList );
+    }
+    CSLDestroy(papszMD_RPC);
+
+    CSLDestroy(papszMD_GEOLOCATION);
+
+    delete poDriver;
+}
+
+/************************************************************************/
+/*                      CloseDependentDatasets()                        */
+/************************************************************************/
+
+int GDALOverviewDataset::CloseDependentDatasets()
+{
+    bool bRet = false;
+
+    if( poMainDS )
+    {
+        for( int i = 0; i < nBands; ++i )
+        {
+            GDALOverviewBand* const band =
+                dynamic_cast<GDALOverviewBand*>(papoBands[i]);
+            if( band == nullptr )
+            {
+                CPLError( CE_Fatal, CPLE_AppDefined,
+                            "OverviewBand cast fail." );
+                return false;
+            }
+            band->poUnderlyingBand = nullptr;
+        }
+        if( poMainDS->ReleaseRef() )
+            bRet = true;
+        poMainDS = nullptr;
+    }
+
+    if( m_poMaskBand )
+    {
+        m_poMaskBand->poUnderlyingBand = nullptr;
+        delete m_poMaskBand;
+        m_poMaskBand = nullptr;
+    }
+
+    return bRet;
+}
+
+/************************************************************************/
+/*                             IRasterIO()                              */
+/*                                                                      */
+/*      The default implementation of IRasterIO() is to pass the        */
+/*      request off to each band objects rasterio methods with          */
+/*      appropriate arguments.                                          */
+/************************************************************************/
+
+CPLErr GDALOverviewDataset::IRasterIO( GDALRWFlag eRWFlag,
+                                       int nXOff, int nYOff,
+                                       int nXSize, int nYSize,
+                                       void * pData,
+                                       int nBufXSize, int nBufYSize,
+                                       GDALDataType eBufType,
+                                       int nBandCount, int *panBandMap,
+                                       GSpacing nPixelSpace,
+                                       GSpacing nLineSpace,
+                                       GSpacing nBandSpace,
+                                       GDALRasterIOExtraArg* psExtraArg )
+
+{
+    // In case the overview bands are really linked to a dataset, then issue
+    // the request to that dataset.
+    if( poOvrDS != nullptr )
+    {
+        return poOvrDS->RasterIO(
+            eRWFlag, nXOff, nYOff, nXSize, nYSize, pData, nBufXSize, nBufYSize,
+            eBufType, nBandCount, panBandMap, nPixelSpace,
+            nLineSpace, nBandSpace,
+            psExtraArg);
+    }
+
+    GDALProgressFunc  pfnProgressGlobal = psExtraArg->pfnProgress;
+    void *pProgressDataGlobal = psExtraArg->pProgressData;
+    CPLErr eErr = CE_None;
+
+    for( int iBandIndex = 0;
+         iBandIndex < nBandCount && eErr == CE_None;
+         ++iBandIndex )
+    {
+        GDALOverviewBand *poBand =
+            dynamic_cast<GDALOverviewBand *>(
+                GetRasterBand(panBandMap[iBandIndex]) );
+        if( poBand == nullptr )
+        {
+            eErr = CE_Failure;
+            break;
+        }
+
+        GByte *pabyBandData =
+            static_cast<GByte *>(pData) + iBandIndex * nBandSpace;
+
+        psExtraArg->pfnProgress = GDALScaledProgress;
+        psExtraArg->pProgressData =
+            GDALCreateScaledProgress( 1.0 * iBandIndex / nBandCount,
+                                      1.0 * (iBandIndex + 1) / nBandCount,
+                                      pfnProgressGlobal,
+                                      pProgressDataGlobal );
+
+        eErr = poBand->IRasterIO( eRWFlag, nXOff, nYOff, nXSize, nYSize,
+                                  pabyBandData,
+                                  nBufXSize, nBufYSize,
+                                  eBufType, nPixelSpace,
+                                  nLineSpace, psExtraArg );
+
+        GDALDestroyScaledProgress( psExtraArg->pProgressData );
+    }
+
+    psExtraArg->pfnProgress = pfnProgressGlobal;
+    psExtraArg->pProgressData = pProgressDataGlobal;
+
+    return eErr;
+}
+
+/************************************************************************/
+/*                           GetSpatialRef()                            */
+/************************************************************************/
+
+const OGRSpatialReference *GDALOverviewDataset::GetSpatialRef() const
+
+{
+    return poMainDS->GetSpatialRef();
+}
+
+/************************************************************************/
+/*                          GetGeoTransform()                           */
+/************************************************************************/
+
+CPLErr GDALOverviewDataset::GetGeoTransform( double * padfTransform )
+
+{
+    double adfGeoTransform[6] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
+    if( poMainDS->GetGeoTransform(adfGeoTransform) != CE_None )
+        return CE_Failure;
+
+    adfGeoTransform[1] *=
+        static_cast<double>(poMainDS->GetRasterXSize()) / nRasterXSize;
+    adfGeoTransform[2] *=
+        static_cast<double>(poMainDS->GetRasterYSize()) / nRasterYSize;
+    adfGeoTransform[4] *=
+        static_cast<double>(poMainDS->GetRasterXSize()) / nRasterXSize;
+    adfGeoTransform[5] *=
+        static_cast<double>(poMainDS->GetRasterYSize()) / nRasterYSize;
+
+    memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
+
+    return CE_None;
+}
+
+/************************************************************************/
+/*                            GetGCPCount()                             */
+/************************************************************************/
+
+int GDALOverviewDataset::GetGCPCount()
+
+{
+    return poMainDS->GetGCPCount();
+}
+
+/************************************************************************/
+/*                          GetGCPSpatialRef()                          */
+/************************************************************************/
+
+const OGRSpatialReference *GDALOverviewDataset::GetGCPSpatialRef() const
+
+{
+    return poMainDS->GetGCPSpatialRef();
+}
+
+/************************************************************************/
+/*                               GetGCPs()                              */
+/************************************************************************/
+
+const GDAL_GCP *GDALOverviewDataset::GetGCPs()
+
+{
+    if( pasGCPList != nullptr )
+        return pasGCPList;
+
+    const GDAL_GCP* pasGCPsMain = poMainDS->GetGCPs();
+    if( pasGCPsMain == nullptr )
+        return nullptr;
+    nGCPCount = poMainDS->GetGCPCount();
+
+    pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPsMain );
+    for( int i = 0; i < nGCPCount; ++i )
+    {
+        pasGCPList[i].dfGCPPixel *= static_cast<double>(nRasterXSize) /
+            poMainDS->GetRasterXSize();
+        pasGCPList[i].dfGCPLine *= static_cast<double>(nRasterYSize) /
+            poMainDS->GetRasterYSize();
+    }
+    return pasGCPList;
+}
+
+/************************************************************************/
+/*                             Rescale()                                */
+/************************************************************************/
+
+void GDALOverviewDataset::Rescale( char**& papszMD, const char* pszItem,
+                                   double dfRatio, double dfDefaultVal )
+{
+    double dfVal =
+        CPLAtofM( CSLFetchNameValueDef(papszMD, pszItem,
+                                       CPLSPrintf("%.18g", dfDefaultVal)) );
+    dfVal *= dfRatio;
+    papszMD = CSLSetNameValue(papszMD, pszItem, CPLSPrintf("%.18g", dfVal));
+}
+
+/************************************************************************/
+/*                            GetMetadata()                             */
+/************************************************************************/
+
+char  **GDALOverviewDataset::GetMetadata( const char * pszDomain )
+{
+    if( poOvrDS != nullptr )
+    {
+        char** papszMD = poOvrDS->GetMetadata(pszDomain);
+        if( papszMD != nullptr )
+            return papszMD;
+    }
+
+    char** papszMD = poMainDS->GetMetadata(pszDomain);
+
+    // We may need to rescale some values from the RPC metadata domain.
+    if( pszDomain != nullptr && EQUAL(pszDomain, MD_DOMAIN_RPC) &&
+        papszMD != nullptr )
+    {
+        if( papszMD_RPC )
+            return papszMD_RPC;
+        papszMD_RPC = CSLDuplicate(papszMD);
+
+        Rescale( papszMD_RPC, RPC_LINE_OFF,
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 0.0 );
+        Rescale( papszMD_RPC, RPC_LINE_SCALE,
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 1.0 );
+        Rescale( papszMD_RPC, RPC_SAMP_OFF,
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 0.0 );
+        Rescale( papszMD_RPC, RPC_SAMP_SCALE,
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 1.0 );
+
+        papszMD = papszMD_RPC;
+    }
+
+    // We may need to rescale some values from the GEOLOCATION metadata domain.
+    if( pszDomain != nullptr && EQUAL(pszDomain, "GEOLOCATION") &&
+        papszMD != nullptr )
+    {
+        if( papszMD_GEOLOCATION )
+            return papszMD_GEOLOCATION;
+        papszMD_GEOLOCATION = CSLDuplicate(papszMD);
+
+        Rescale( papszMD_GEOLOCATION, "PIXEL_OFFSET",
+                 static_cast<double>(poMainDS->GetRasterXSize()) /
+                 nRasterXSize, 0.0 );
+        Rescale( papszMD_GEOLOCATION, "LINE_OFFSET",
+                 static_cast<double>(poMainDS->GetRasterYSize()) /
+                 nRasterYSize, 0.0 );
+
+        Rescale( papszMD_GEOLOCATION, "PIXEL_STEP",
+                 static_cast<double>(nRasterXSize) / poMainDS->GetRasterXSize(),
+                 1.0 );
+        Rescale( papszMD_GEOLOCATION, "LINE_STEP",
+                 static_cast<double>(nRasterYSize) / poMainDS->GetRasterYSize(),
+                 1.0 );
+
+        papszMD = papszMD_GEOLOCATION;
+    }
+
+    return papszMD;
+}
+
+/************************************************************************/
+/*                          GetMetadataItem()                           */
+/************************************************************************/
+
+const char *GDALOverviewDataset::GetMetadataItem( const char * pszName,
+                                                  const char * pszDomain )
+{
+    if( poOvrDS != nullptr )
+    {
+        const char* pszValue = poOvrDS->GetMetadataItem(pszName, pszDomain);
+        if( pszValue != nullptr )
+            return pszValue;
+    }
+
+    if( pszDomain != nullptr && (EQUAL(pszDomain, "RPC") ||
+                              EQUAL(pszDomain, "GEOLOCATION")) )
+    {
+        char** papszMD = GetMetadata(pszDomain);
+        return CSLFetchNameValue(papszMD, pszName);
+    }
+
+    return poMainDS->GetMetadataItem(pszName, pszDomain);
+}
+
+/************************************************************************/
+/*                          GDALOverviewBand()                          */
+/************************************************************************/
+
+GDALOverviewBand::GDALOverviewBand( GDALOverviewDataset* poDSIn, int nBandIn )
+{
+    poDS = poDSIn;
+    nBand = nBandIn;
+    nRasterXSize = poDSIn->nRasterXSize;
+    nRasterYSize = poDSIn->nRasterYSize;
+    if( nBandIn == 0 )
+    {
+        poUnderlyingBand = poDSIn->poMainDS->GetRasterBand(1)->
+                            GetOverview(poDSIn->nOvrLevel)->GetMaskBand();
+    }
+    else
+    {
+        poUnderlyingBand = poDSIn->poMainDS->GetRasterBand(nBandIn)->
+                            GetOverview(poDSIn->nOvrLevel);
+    }
+    eDataType = poUnderlyingBand->GetRasterDataType();
+    poUnderlyingBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
+}
+
+/************************************************************************/
+/*                         ~GDALOverviewBand()                          */
+/************************************************************************/
+
+GDALOverviewBand::~GDALOverviewBand()
+{
+    GDALOverviewBand::FlushCache();
+}
+
+/************************************************************************/
+/*                              FlushCache()                            */
+/************************************************************************/
+
+CPLErr GDALOverviewBand::FlushCache()
+{
+    if( poUnderlyingBand )
+        return poUnderlyingBand->FlushCache();
+    return CE_None;
+}
+
+/************************************************************************/
+/*                        RefUnderlyingRasterBand()                     */
+/************************************************************************/
+
+GDALRasterBand* GDALOverviewBand::RefUnderlyingRasterBand()
+{
+    if( poUnderlyingBand )
+        return poUnderlyingBand;
+
+    return nullptr;
+}
+
+/************************************************************************/
+/*                         GetOverviewCount()                           */
+/************************************************************************/
+
+int GDALOverviewBand::GetOverviewCount()
+{
+    GDALOverviewDataset * const poOvrDS =
+        dynamic_cast<GDALOverviewDataset *>(poDS);
+    if( poOvrDS == nullptr )
+    {
+        CPLError( CE_Fatal, CPLE_AppDefined, "OverviewDataset cast fail." );
+        return 0;
+    }
+    if( poOvrDS->bThisLevelOnly )
+        return 0;
+    GDALDataset * const poMainDS = poOvrDS->poMainDS;
+    GDALRasterBand* poMainBand =
+        ( nBand == 0 ) ? poMainDS->GetRasterBand(1)->GetMaskBand() :
+                         poMainDS->GetRasterBand(nBand);
+    return poMainBand->GetOverviewCount() - poOvrDS->nOvrLevel - 1;
+}
+
+/************************************************************************/
+/*                           GetOverview()                              */
+/************************************************************************/
+
+GDALRasterBand *GDALOverviewBand::GetOverview( int iOvr )
+{
+    if( iOvr < 0 || iOvr >= GetOverviewCount() )
+        return nullptr;
+    GDALOverviewDataset * const poOvrDS =
+        dynamic_cast<GDALOverviewDataset *>(poDS);
+    if( poOvrDS == nullptr )
+    {
+        CPLError( CE_Fatal, CPLE_AppDefined, "OverviewDataset cast fail." );
+        return nullptr;
+    }
+    GDALDataset * const poMainDS = poOvrDS->poMainDS;
+    GDALRasterBand* poMainBand =
+        ( nBand == 0 ) ? poMainDS->GetRasterBand(1)->GetMaskBand() :
+                         poMainDS->GetRasterBand(nBand);
+    return poMainBand->GetOverview(iOvr + poOvrDS->nOvrLevel + 1);
+}
+
+/************************************************************************/
+/*                           GetMaskFlags()                             */
+/************************************************************************/
+
+int GDALOverviewBand::GetMaskFlags()
+{
+    GDALOverviewDataset * const poOvrDS = cpl::down_cast<GDALOverviewDataset *>(poDS);
+    if( nBand != 0 && poOvrDS->m_poMaskBand )
+        return GMF_PER_DATASET;
+    return GDALProxyRasterBand::GetMaskFlags();
+}
+
+/************************************************************************/
+/*                           GetMaskBand()                              */
+/************************************************************************/
+
+GDALRasterBand* GDALOverviewBand::GetMaskBand()
+{
+    GDALOverviewDataset * const poOvrDS = cpl::down_cast<GDALOverviewDataset *>(poDS);
+    if( nBand != 0 && poOvrDS->m_poMaskBand )
+        return poOvrDS->m_poMaskBand;
+    return GDALProxyRasterBand::GetMaskBand();
+}
diff --git a/src/types.hpp b/src/types.hpp
index 9fc93f3..7c61ad2 100644
--- a/src/types.hpp
+++ b/src/types.hpp
@@ -25,6 +25,7 @@
 #include <cstdint>              // uint16_t
 
 #include "Bounds.hpp"
+#include "Coordinate3D.hpp"
 
 /// All terrain related data types reside in this namespace
 namespace ctb {
@@ -38,7 +39,8 @@ namespace ctb {
   // Complex types
   typedef Bounds<i_tile> TileBounds;      ///< Tile extents in tile coordinates
   typedef Coordinate<i_pixel> PixelPoint; ///< The location of a pixel
-  typedef Coordinate<double> CRSPoint; ///< A Coordinate Reference System coordinate
+  typedef Coordinate<double> CRSPoint;    ///< A Coordinate Reference System coordinate
+  typedef Coordinate3D<double> CRSVertex; ///< A 3D-Vertex of a mesh or tile in CRS coordinates
   typedef Bounds<double> CRSBounds;       ///< Extents in CRS coordinates
   typedef Coordinate<i_tile> TilePoint;   ///< The location of a tile
 
diff --git a/tools/ctb-tile.cpp b/tools/ctb-tile.cpp
index 3b4b679..1bce7ad 100644
--- a/tools/ctb-tile.cpp
+++ b/tools/ctb-tile.cpp
@@ -50,6 +50,9 @@
 #include "GlobalMercator.hpp"
 #include "RasterIterator.hpp"
 #include "TerrainIterator.hpp"
+#include "MeshIterator.hpp"
+#include "GDALDatasetReader.hpp"
+#include "CTBFileTileSerializer.hpp"
 
 using namespace std;
 using namespace ctb;
@@ -73,7 +76,11 @@ class TerrainBuild : public Command {
     startZoom(-1),
     endZoom(-1),
     verbosity(1),
-    resume(false)
+    resume(false),
+    meshQualityFactor(1.0),
+    metadata(false),
+    cesiumFriendly(false),
+    vertexNormals(false)
   {}
 
   void
@@ -198,6 +205,26 @@ class TerrainBuild : public Command {
     return  (command->argc == 1) ? command->argv[0] : NULL;
   }
 
+  static void
+    setMeshQualityFactor(command_t *command) {
+    static_cast<TerrainBuild *>(Command::self(command))->meshQualityFactor = atof(command->arg);
+  }
+
+  static void
+    setMetadata(command_t *command) {
+    static_cast<TerrainBuild *>(Command::self(command))->metadata = true;
+  }
+
+  static void
+    setCesiumFriendly(command_t *command) {
+    static_cast<TerrainBuild *>(Command::self(command))->cesiumFriendly = true;
+  }
+
+  static void
+    setVertexNormals(command_t *command) {
+    static_cast<TerrainBuild *>(Command::self(command))->vertexNormals = true;
+  }
+
   const char *outputDir,
     *outputFormat,
     *profile;
@@ -212,53 +239,12 @@ class TerrainBuild : public Command {
 
   CPLStringList creationOptions;
   TilerOptions tilerOptions;
-};
-
-/**
- * Create a filename for a tile coordinate
- *
- * This also creates the tile directory structure.
- */
-static string
-getTileFilename(const TileCoordinate *coord, const string dirname, const char *extension) {
-  static mutex mutex;
-  VSIStatBufL stat;
-  string filename = concat(dirname, coord->zoom, osDirSep, coord->x);
 
-  lock_guard<std::mutex> lock(mutex);
-
-  // Check whether the `{zoom}/{x}` directory exists or not
-  if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) {
-    filename = concat(dirname, coord->zoom);
-
-    // Check whether the `{zoom}` directory exists or not
-    if (VSIStatExL(filename.c_str(), &stat, VSI_STAT_EXISTS_FLAG | VSI_STAT_NATURE_FLAG)) {
-      // Create the `{zoom}` directory
-      if (VSIMkdir(filename.c_str(), 0755))
-        throw CTBException("Could not create the zoom level directory");
-
-    } else if (!VSI_ISDIR(stat.st_mode)) {
-      throw CTBException("Zoom level file path is not a directory");
-    }
-
-    // Create the `{zoom}/{x}` directory
-    filename += concat(osDirSep, coord->x);
-    if (VSIMkdir(filename.c_str(), 0755))
-      throw CTBException("Could not create the x level directory");
-
-  } else if (!VSI_ISDIR(stat.st_mode)) {
-    throw CTBException("X level file path is not a directory");
-  }
-
-  // Create the filename itself, adding the extension if required
-  filename += concat(osDirSep, coord->y);
-  if (extension != NULL) {
-    filename += ".";
-    filename += extension;
-  }
-
-  return filename;
-}
+  double meshQualityFactor;
+  bool metadata;
+  bool cesiumFriendly;
+  bool vertexNormals;
+};
 
 /**
  * Increment a TilerIterator whilst cooperating between threads
@@ -269,9 +255,9 @@ getTileFilename(const TileCoordinate *coord, const string dirname, const char *e
  * to ensure all tiles are iterated over consecutively.  It assumes individual
  * tile iterators point to the same source GDAL dataset.
  */
+static int globalIteratorIndex = 0; // keep track of where we are globally
 template<typename T> int
 incrementIterator(T &iter, int currentIndex) {
-  static int globalIteratorIndex = 0; // keep track of where we are globally
   static mutex mutex;        // ensure iterations occur serially between threads
 
   lock_guard<std::mutex> lock(mutex);
@@ -332,6 +318,10 @@ showProgress(int currentIndex, string filename) {
 
   return progressFunc(currentIndex / (double) iteratorSize, message.c_str(), NULL);
 }
+int
+showProgress(int currentIndex) {
+  return progressFunc(currentIndex / (double) iteratorSize, NULL, NULL);
+}
 
 static bool
 fileExists(const std::string& filename) {
@@ -339,9 +329,250 @@ fileExists(const std::string& filename) {
   return VSIStatExL(filename.c_str(), &statbuf, VSI_STAT_EXISTS_FLAG) == 0;
 }
 
+static bool
+fileCopy(const std::string& sourceFilename, const std::string& targetFilename) {
+  FILE *fp_s = VSIFOpen(sourceFilename.c_str(), "rb");
+  if (!fp_s) return false;
+  FILE *fp_t = VSIFOpen(targetFilename.c_str(), "wb");
+  if (!fp_t) return false;
+
+  VSIFSeek(fp_s, 0, SEEK_END);
+  long fileSize = VSIFTell(fp_s);
+
+  if (fileSize > 0)
+  {
+    VSIFSeek(fp_s, 0, SEEK_SET);
+
+    void* buffer = VSIMalloc(fileSize);
+    VSIFRead(buffer, 1, fileSize, fp_s);
+    VSIFWrite(buffer, 1, fileSize, fp_t);
+    VSIFree(buffer);
+  }
+  VSIFClose(fp_t);
+  VSIFClose(fp_s);
+
+  return fileSize > 0;
+}
+
+/// Handle the terrain metadata
+class TerrainMetadata {
+public:
+  TerrainMetadata() {
+  }
+
+  // Defines the valid tile indexes of a level in a Tileset
+  struct LevelInfo {
+  public:
+    LevelInfo() {
+      startX = startY = std::numeric_limits<int>::max();
+      finalX = finalY = std::numeric_limits<int>::min();
+    }
+    int startX, startY;
+    int finalX, finalY;
+
+    inline void add(const TileCoordinate *coordinate) {
+      startX = std::min(startX, (int)coordinate->x);
+      startY = std::min(startY, (int)coordinate->y);
+      finalX = std::max(finalX, (int)coordinate->x);
+      finalY = std::max(finalY, (int)coordinate->y);
+    }
+    inline void add(const LevelInfo &level) {
+      startX = std::min(startX, level.startX);
+      startY = std::min(startY, level.startY);
+      finalX = std::max(finalX, level.finalX);
+      finalY = std::max(finalY, level.finalY);
+    }
+  };
+  std::vector<LevelInfo> levels;
+
+  // Defines the bounding box covered by the Terrain
+  CRSBounds bounds;
+
+  // Add metadata of the specified Coordinate
+  void add(const Grid &grid, const TileCoordinate *coordinate) {
+    CRSBounds tileBounds = grid.tileBounds(*coordinate);
+    i_zoom zoom = coordinate->zoom;
+
+    if ((1 + zoom) > levels.size()) {
+      levels.resize(1 + zoom, LevelInfo());
+    }
+    LevelInfo &level = levels[zoom];
+    level.add(coordinate);
+
+    if (bounds.getMaxX() == bounds.getMinX()) {
+      bounds = tileBounds;
+    }
+    else {
+      bounds.setMinX(std::min(bounds.getMinX(), tileBounds.getMinX()));
+      bounds.setMinY(std::min(bounds.getMinY(), tileBounds.getMinY()));
+      bounds.setMaxX(std::max(bounds.getMaxX(), tileBounds.getMaxX()));
+      bounds.setMaxY(std::max(bounds.getMaxY(), tileBounds.getMaxY()));
+    }
+  }
+
+  // Add metadata info
+  void add(const TerrainMetadata &otherMetadata) {
+    if (otherMetadata.levels.size() > 0) {
+      const CRSBounds &otherBounds = otherMetadata.bounds;
+
+      if(otherMetadata.levels.size() > levels.size())
+      {
+        levels.resize(otherMetadata.levels.size(), LevelInfo());
+      }
+      for (size_t i = 0; i < otherMetadata.levels.size(); i++) {
+        levels[i].add(otherMetadata.levels[i]);
+      }
+
+      bounds.setMinX(std::min(bounds.getMinX(), otherBounds.getMinX()));
+      bounds.setMinY(std::min(bounds.getMinY(), otherBounds.getMinY()));
+      bounds.setMaxX(std::max(bounds.getMaxX(), otherBounds.getMaxX()));
+      bounds.setMaxY(std::max(bounds.getMaxY(), otherBounds.getMaxY()));
+    }
+  }
+
+  /// Output the layer.json metadata file
+  /// http://help.agi.com/TerrainServer/RESTAPIGuide.html
+  /// Example:
+  /// https://assets.agi.com/stk-terrain/v1/tilesets/world/tiles/layer.json
+  void writeJsonFile(const std::string &filename, const std::string &datasetName, const std::string &outputFormat = "Terrain", const std::string &profile = "geodetic", bool writeVertexNormals = false) const {
+    FILE *fp = fopen(filename.c_str(), "w");
+
+    if (fp == NULL) {
+      throw CTBException("Failed to open metadata file");
+    }
+
+    fprintf(fp, "{\n");
+    fprintf(fp, "  \"tilejson\": \"2.1.0\",\n");
+    fprintf(fp, "  \"name\": \"%s\",\n", datasetName.c_str());
+    fprintf(fp, "  \"description\": \"\",\n");
+    fprintf(fp, "  \"version\": \"1.1.0\",\n");
+
+    if (strcmp(outputFormat.c_str(), "Terrain") == 0) {
+      fprintf(fp, "  \"format\": \"heightmap-1.0\",\n");
+    }
+    else if (strcmp(outputFormat.c_str(), "Mesh") == 0) {
+      fprintf(fp, "  \"format\": \"quantized-mesh-1.0\",\n");
+    }
+    else {
+      fprintf(fp, "  \"format\": \"GDAL\",\n");
+    }
+    fprintf(fp, "  \"attribution\": \"\",\n");
+    fprintf(fp, "  \"schema\": \"tms\",\n");
+    if (writeVertexNormals) {
+      fprintf(fp, "  \"extensions\": [ \"octvertexnormals\" ],\n");
+    }
+    fprintf(fp, "  \"tiles\": [ \"{z}/{x}/{y}.terrain?v={version}\" ],\n");
+
+    if (strcmp(profile.c_str(), "geodetic") == 0) {
+      fprintf(fp, "  \"projection\": \"EPSG:4326\",\n");
+    }
+    else {
+      fprintf(fp, "  \"projection\": \"EPSG:3857\",\n");
+    }
+    fprintf(fp, "  \"bounds\": [ %.2f, %.2f, %.2f, %.2f ],\n",
+      bounds.getMinX(),
+      bounds.getMinY(),
+      bounds.getMaxX(),
+      bounds.getMaxY());
+
+    fprintf(fp, "  \"available\": [\n");
+    for (size_t i = 0, icount = levels.size(); i < icount; i++) {
+      const LevelInfo &level = levels[i];
+
+      if (i > 0)
+        fprintf(fp, "   ,[ ");
+      else
+        fprintf(fp, "    [ ");
+
+      if (level.finalX >= level.startX) {
+        fprintf(fp, "{ \"startX\": %li, \"startY\": %li, \"endX\": %li, \"endY\": %li }",
+          level.startX,
+          level.startY,
+          level.finalX,
+          level.finalY);
+      }
+      fprintf(fp, " ]\n");
+    }
+    fprintf(fp, "  ]\n");
+
+    fprintf(fp, "}\n");
+    fclose(fp);
+  }
+};
+
+/// Create an empty root temporary elevation file (GTiff)
+static std::string 
+createEmptyRootElevationFile(std::string &fileName, const Grid &grid, const TileCoordinate& coord) {
+  GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");
+
+  if (poDriver == NULL) {
+    throw CTBException("Could not retrieve GTiff GDAL driver");
+  }
+
+  // Create the geo transform for this temporary elevation tile.
+  // We apply an 1-degree negative offset to avoid problems in borders.
+  CRSBounds tileBounds = grid.tileBounds(coord);
+  tileBounds.setMinX(tileBounds.getMinX() + 1);
+  tileBounds.setMinY(tileBounds.getMinY() + 1);
+  tileBounds.setMaxX(tileBounds.getMaxX() - 1);
+  tileBounds.setMaxY(tileBounds.getMaxY() - 1);
+  const i_tile tileSize = grid.tileSize() - 2;
+  const double resolution = tileBounds.getWidth() / tileSize;
+  double adfGeoTransform[6] = { tileBounds.getMinX(), resolution, 0, tileBounds.getMaxY(), 0, -resolution };
+
+  // Create the spatial reference system for the file
+  OGRSpatialReference oSRS;
+
+  #if ( GDAL_VERSION_MAJOR >= 3 )
+  oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
+  #endif
+  
+  if (oSRS.importFromEPSG(4326) != OGRERR_NONE) {
+    throw CTBException("Could not create EPSG:4326 spatial reference");
+  }
+  char *pszDstWKT = NULL;
+  if (oSRS.exportToWkt(&pszDstWKT) != OGRERR_NONE) {
+    CPLFree(pszDstWKT);
+    throw CTBException("Could not create EPSG:4326 WKT string");
+  }
+
+  // Create the GTiff file
+  fileName += ".tif";
+  GDALDataset *poDataset = poDriver->Create(fileName.c_str(), tileSize, tileSize, 1, GDT_Float32, NULL);
+
+  // Set the projection
+  if (poDataset->SetProjection(pszDstWKT) != CE_None) {
+    poDataset->Release();
+    CPLFree(pszDstWKT);
+    throw CTBException("Could not set projection on temporary elevation file");
+  }
+  CPLFree(pszDstWKT);
+
+  // Apply the geo transform
+  if (poDataset->SetGeoTransform(adfGeoTransform) != CE_None) {
+    poDataset->Release();
+    throw CTBException("Could not set projection on temporary elevation file");
+  }
+
+  // Finally write the height data
+  float *rasterHeights = (float *)CPLCalloc(tileSize * tileSize, sizeof(float));
+  GDALRasterBand *heightsBand = poDataset->GetRasterBand(1);
+  if (heightsBand->RasterIO(GF_Write, 0, 0, tileSize, tileSize, 
+                            (void *)rasterHeights, tileSize, tileSize, GDT_Float32, 
+                            0, 0) != CE_None) {
+    CPLFree(rasterHeights);
+    throw CTBException("Could not write heights on temporary elevation file");
+  }
+  CPLFree(rasterHeights);
+
+  poDataset->FlushCache();
+  poDataset->Release();
+  return fileName;
+}
+
 /// Output GDAL tiles represented by a tiler to a directory
 static void
-buildGDAL(const RasterTiler &tiler, TerrainBuild *command) {
+buildGDAL(GDALSerializer &serializer, const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) {
   GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName(command->outputFormat);
 
   if (poDriver == NULL) {
@@ -353,7 +584,6 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command) {
   }
 
   const char *extension = poDriver->GetMetadataItem(GDAL_DMD_EXTENSION);
-  const string dirname = string(command->outputDir) + osDirSep;
   i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom,
     endZoom = (command->endZoom < 0) ? 0 : command->endZoom;
 
@@ -363,61 +593,109 @@ buildGDAL(const RasterTiler &tiler, TerrainBuild *command) {
 
   while (!iter.exhausted()) {
     const TileCoordinate *coordinate = iter.GridIterator::operator*();
-    GDALDataset *poDstDS;
-    const string filename = getTileFilename(coordinate, dirname, extension);
+    if (metadata) metadata->add(tiler.grid(), coordinate);
 
-    
-    if( !command->resume || !fileExists(filename) ) {
+    if (serializer.mustSerializeCoordinate(coordinate)) {
       GDALTile *tile = *iter;
-      const string temp_filename = concat(filename, ".tmp");
-      poDstDS = poDriver->CreateCopy(temp_filename.c_str(), tile->dataset, FALSE,
-                                     command->creationOptions.List(), NULL, NULL );
+      serializer.serializeTile(tile, poDriver, extension, command->creationOptions);
       delete tile;
-
-      // Close the datasets, flushing data to destination
-      if (poDstDS == NULL) {
-        throw CTBException("Could not create GDAL tile");
-      }
-
-      GDALClose(poDstDS);
-
-      if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) {
-        throw new CTBException("Could not rename temporary file");
-      }
     }
 
     currentIndex = incrementIterator(iter, currentIndex);
-    showProgress(currentIndex, filename);
+    showProgress(currentIndex);
   }
 }
 
 /// Output terrain tiles represented by a tiler to a directory
 static void
-buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) {
-  const string dirname = string(command->outputDir) + osDirSep;
+buildTerrain(TerrainSerializer &serializer, const TerrainTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) {
   i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom,
     endZoom = (command->endZoom < 0) ? 0 : command->endZoom;
 
   TerrainIterator iter(tiler, startZoom, endZoom);
   int currentIndex = incrementIterator(iter, 0);
   setIteratorSize(iter);
+  GDALDatasetReaderWithOverviews reader(tiler);
 
   while (!iter.exhausted()) {
     const TileCoordinate *coordinate = iter.GridIterator::operator*();
-    const string filename = getTileFilename(coordinate, dirname, "terrain");
-
-    if( !command->resume || !fileExists(filename) ) {
-      TerrainTile *tile = *iter;
-      const string temp_filename = concat(filename, ".tmp");
+    if (metadata) metadata->add(tiler.grid(), coordinate);
 
-      tile->writeFile(temp_filename.c_str());
+    if (serializer.mustSerializeCoordinate(coordinate)) {
+      TerrainTile *tile = iter.operator*(&reader);
+      serializer.serializeTile(tile);
       delete tile;
+    }
 
-      if (VSIRename(temp_filename.c_str(), filename.c_str()) != 0) {
-        throw new CTBException("Could not rename temporary file");
-      }
+    currentIndex = incrementIterator(iter, currentIndex);
+    showProgress(currentIndex);
+  }
+}
+
+/// Output mesh tiles represented by a tiler to a directory
+static void
+buildMesh(MeshSerializer &serializer, const MeshTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata, bool writeVertexNormals = false) {
+  i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom,
+    endZoom = (command->endZoom < 0) ? 0 : command->endZoom;
+
+  // DEBUG Chunker:
+  #if 0
+  const string dirname = string(command->outputDir) + osDirSep;
+  TileCoordinate coordinate(13, 8102, 6047);
+  MeshTile *tile = tiler.createMesh(tiler.dataset(), coordinate);
+  //
+  const string txtname = CTBFileTileSerializer::getTileFilename(&coordinate, dirname, "wkt");
+  const Mesh &mesh = tile->getMesh();
+  mesh.writeWktFile(txtname.c_str());
+  //
+  CRSBounds bounds = tiler.grid().tileBounds(coordinate);
+  double x = bounds.getMinX() + 0.5 * (bounds.getMaxX() - bounds.getMinX());
+  double y = bounds.getMinY() + 0.5 * (bounds.getMaxY() - bounds.getMinY());
+  CRSPoint point(x,y);
+  TileCoordinate c = tiler.grid().crsToTile(point, coordinate.zoom);
+  //
+  const string filename = CTBFileTileSerializer::getTileFilename(&coordinate, dirname, "terrain");
+  tile->writeFile(filename.c_str(), writeVertexNormals);
+  delete tile;
+  return;
+  #endif
+
+  MeshIterator iter(tiler, startZoom, endZoom);
+  int currentIndex = incrementIterator(iter, 0);
+  setIteratorSize(iter);
+  GDALDatasetReaderWithOverviews reader(tiler);
+
+  while (!iter.exhausted()) {
+    const TileCoordinate *coordinate = iter.GridIterator::operator*();
+    if (metadata) metadata->add(tiler.grid(), coordinate);
+
+    if (serializer.mustSerializeCoordinate(coordinate)) {
+      MeshTile *tile = iter.operator*(&reader);
+      serializer.serializeTile(tile, writeVertexNormals);
+      delete tile;
     }
 
+    currentIndex = incrementIterator(iter, currentIndex);
+    showProgress(currentIndex);
+  }
+}
+
+static void
+buildMetadata(const RasterTiler &tiler, TerrainBuild *command, TerrainMetadata *metadata) {
+  const string dirname = string(command->outputDir) + osDirSep;
+  i_zoom startZoom = (command->startZoom < 0) ? tiler.maxZoomLevel() : command->startZoom,
+    endZoom = (command->endZoom < 0) ? 0 : command->endZoom;
+
+  const std::string filename = concat(dirname, "layer.json"); 
+
+  RasterIterator iter(tiler, startZoom, endZoom);
+  int currentIndex = incrementIterator(iter, 0);
+  setIteratorSize(iter);
+
+  while (!iter.exhausted()) {
+    const TileCoordinate *coordinate = iter.GridIterator::operator*();
+    if (metadata) metadata->add(tiler.grid(), coordinate);
+
     currentIndex = incrementIterator(iter, currentIndex);
     showProgress(currentIndex, filename);
   }
@@ -429,28 +707,51 @@ buildTerrain(const TerrainTiler &tiler, TerrainBuild *command) {
  * This function is designed to be run in a separate thread.
  */
 static int
-runTiler(TerrainBuild *command, Grid *grid) {
-  GDALDataset  *poDataset = (GDALDataset *) GDALOpen(command->getInputFilename(), GA_ReadOnly);
+runTiler(const char *inputFilename, TerrainBuild *command, Grid *grid, TerrainMetadata *metadata) {
+  GDALDataset  *poDataset = (GDALDataset *) GDALOpen(inputFilename, GA_ReadOnly);
   if (poDataset == NULL) {
     cerr << "Error: could not open GDAL dataset" << endl;
     return 1;
   }
 
+  // Metadata of only this thread, it will be joined to global later
+  TerrainMetadata *threadMetadata = metadata ? new TerrainMetadata() : NULL;
+
+  // Choose serializer of tiles (Directory of files, MBTiles store...)
+  CTBFileTileSerializer serializer(string(command->outputDir) + osDirSep, command->resume);
+
   try {
-    if (strcmp(command->outputFormat, "Terrain") == 0) {
+    serializer.startSerialization();
+
+    if (command->metadata) {
+      const RasterTiler tiler(poDataset, *grid, command->tilerOptions);
+      buildMetadata(tiler, command, threadMetadata);
+    } else if (strcmp(command->outputFormat, "Terrain") == 0) {
       const TerrainTiler tiler(poDataset, *grid);
-      buildTerrain(tiler, command);
+      buildTerrain(serializer, tiler, command, threadMetadata);
+    } else if (strcmp(command->outputFormat, "Mesh") == 0) {
+      const MeshTiler tiler(poDataset, *grid, command->tilerOptions, command->meshQualityFactor);
+      buildMesh(serializer, tiler, command, threadMetadata, command->vertexNormals);
     } else {                    // it's a GDAL format
       const RasterTiler tiler(poDataset, *grid, command->tilerOptions);
-      buildGDAL(tiler, command);
+      buildGDAL(serializer, tiler, command, threadMetadata);
     }
 
   } catch (CTBException &e) {
     cerr << "Error: " << e.what() << endl;
   }
+  serializer.endSerialization();
 
   GDALClose(poDataset);
 
+  // Pass metadata to global instance.
+  if (threadMetadata) {
+    static std::mutex mutex;
+    std::lock_guard<std::mutex> lock(mutex);
+
+    metadata->add(*threadMetadata);
+    delete threadMetadata;
+  }
   return 0;
 }
 
@@ -460,7 +761,7 @@ main(int argc, char *argv[]) {
   TerrainBuild command = TerrainBuild(argv[0], version.cstr);
   command.setUsage("[options] GDAL_DATASOURCE");
   command.option("-o", "--output-dir <dir>", "specify the output directory for the tiles (defaults to working directory)", TerrainBuild::setOutputDir);
-  command.option("-f", "--output-format <format>", "specify the output format for the tiles. This is either `Terrain` (the default) or any format listed by `gdalinfo --formats`", TerrainBuild::setOutputFormat);
+  command.option("-f", "--output-format <format>", "specify the output format for the tiles. This is either `Terrain` (the default), `Mesh` (Chunked LOD mesh), or any format listed by `gdalinfo --formats`", TerrainBuild::setOutputFormat);
   command.option("-p", "--profile <profile>", "specify the TMS profile for the tiles. This is either `geodetic` (the default) or `mercator`", TerrainBuild::setProfile);
   command.option("-c", "--thread-count <count>", "specify the number of threads to use for tile generation. On multicore machines this defaults to the number of CPUs", TerrainBuild::setThreadCount);
   command.option("-t", "--tile-size <size>", "specify the size of the tiles in pixels. This defaults to 65 for terrain tiles and 256 for other GDAL formats", TerrainBuild::setTileSize);
@@ -471,6 +772,10 @@ main(int argc, char *argv[]) {
   command.option("-z", "--error-threshold <threshold>", "specify the error threshold in pixel units for transformation approximation. Larger values should mean faster transforms. Defaults to 0.125", TerrainBuild::setErrorThreshold);
   command.option("-m", "--warp-memory <bytes>", "The memory limit in bytes used for warp operations. Higher settings should be faster. Defaults to a conservative GDAL internal setting.", TerrainBuild::setWarpMemory);
   command.option("-R", "--resume", "Do not overwrite existing files", TerrainBuild::setResume);
+  command.option("-g", "--mesh-qfactor <factor>", "specify the factor to multiply the estimated geometric error to convert heightmaps to irregular meshes. Larger values should mean minor quality. Defaults to 1.0", TerrainBuild::setMeshQualityFactor);
+  command.option("-l", "--layer", "only output the layer.json metadata file", TerrainBuild::setMetadata);
+  command.option("-C", "--cesium-friendly", "Force the creation of missing root tiles to be CesiumJS-friendly", TerrainBuild::setCesiumFriendly);
+  command.option("-N", "--vertex-normals", "Write 'Oct-Encoded Per-Vertex Normals' for Terrain Lighting, only for `Mesh` format", TerrainBuild::setVertexNormals);
   command.option("-q", "--quiet", "only output errors", TerrainBuild::setQuiet);
   command.option("-v", "--verbose", "be more noisy", TerrainBuild::setVerbose);
 
@@ -514,11 +819,16 @@ main(int argc, char *argv[]) {
   vector<future<int>> tasks;
   int threadCount = (command.threadCount > 0) ? command.threadCount : CPLGetNumCPUs();
 
+  // Calculate metadata?
+  const string dirname = string(command.outputDir) + osDirSep;
+  const std::string filename = concat(dirname, "layer.json");
+  TerrainMetadata *metadata = command.metadata ? new TerrainMetadata() : NULL;
+
   // Instantiate the threads using futures from a packaged_task
   for (int i = 0; i < threadCount ; ++i) {
-    packaged_task<int(TerrainBuild *, Grid *)> task(runTiler); // wrap the function
-    tasks.push_back(task.get_future());                        // get a future
-    thread(move(task), &command, &grid).detach(); // launch on a thread
+    packaged_task<int(const char *, TerrainBuild *, Grid *, TerrainMetadata *)> task(runTiler); // wrap the function
+    tasks.push_back(task.get_future()); // get a future
+    thread(move(task), command.getInputFilename(), &command, &grid, metadata).detach(); // launch on a thread
   }
 
   // Synchronise the completion of the threads
@@ -531,8 +841,66 @@ main(int argc, char *argv[]) {
     int retval = task.get();
 
     // return on the first encountered problem
-    if (retval)
+    if (retval) {
+      delete metadata;
       return retval;
+    }
+  }
+
+  // CesiumJS friendly?
+  if (command.cesiumFriendly && (strcmp(command.profile, "geodetic") == 0) && command.endZoom <= 0) {
+
+    // Create missing root tiles if it is necessary
+    if (!command.metadata) {
+      std::string dirName0 = string(command.outputDir) + osDirSep + "0" + osDirSep + "0";
+      std::string dirName1 = string(command.outputDir) + osDirSep + "0" + osDirSep + "1";
+      std::string tileName0 = dirName0 + osDirSep + "0.terrain";
+      std::string tileName1 = dirName1 + osDirSep + "0.terrain";
+
+      i_zoom missingZoom = 65535;
+      ctb::TileCoordinate missingTileCoord(missingZoom, 0, 0);
+      std::string missingTileName;
+
+      if (fileExists(tileName0) && !fileExists(tileName1)) {
+        VSIMkdir(dirName1.c_str(), 0755);
+        missingTileCoord = ctb::TileCoordinate(0, 1, 0);
+        missingTileName = tileName1;
+      }
+      else
+      if (!fileExists(tileName0) && fileExists(tileName1)) {
+        VSIMkdir(dirName0.c_str(), 0755);
+        missingTileCoord = ctb::TileCoordinate(0, 0, 0);
+        missingTileName = tileName0;
+      }
+      if (missingTileCoord.zoom != missingZoom) {
+        globalIteratorIndex = 0; // reset global iterator index
+        command.startZoom = 0;
+        command.endZoom = 0;
+        missingTileName = createEmptyRootElevationFile(missingTileName, grid, missingTileCoord);
+        runTiler (missingTileName.c_str(), &command, &grid, NULL);
+        VSIUnlink(missingTileName.c_str());
+      }
+    }
+
+    // Fix available indexes.
+    if (metadata && metadata->levels.size() > 0) {
+      TerrainMetadata::LevelInfo &level = metadata->levels.at(0);
+      level.startX = 0;
+      level.startY = 0;
+      level.finalX = 1;
+      level.finalY = 0;
+    }
+  }
+
+  // Write Json metadata file?
+  if (metadata) {
+    std::string datasetName(command.getInputFilename());
+    datasetName = datasetName.substr(datasetName.find_last_of("/\\") + 1);
+    const size_t rfindpos = datasetName.rfind('.');
+    if (std::string::npos != rfindpos) datasetName = datasetName.erase(rfindpos);
+
+    metadata->writeJsonFile(filename, datasetName, std::string(command.outputFormat), std::string(command.profile), command.vertexNormals);
+    delete metadata;
   }
 
   return 0;