diff --git a/CMakeLists.txt b/CMakeLists.txt
index 377607f16b..cd6f9fa1dc 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -113,6 +113,16 @@ endif()
 
 find_package(Boost 1.72.0 REQUIRED)
 
+
+find_package(SQLite3)
+if(SQLite3_FOUND)
+    set(NGEN_WITH_SQLITE3 ON)
+    add_compile_definitions(NGEN_WITH_SQLITE3)
+else()
+    set(NGEN_WITH_SQLITE3 OFF)
+endif()
+message(INFO " SQLite3 support is ${NGEN_WITH_SQLITE3}")
+
 # UDUNITS
 # Since UDUNITS is currently not really optional (yet) let's make it default to on...
 if (NOT (DEFINED UDUNITS_ACTIVE))
@@ -259,6 +269,12 @@ endif()
 add_subdirectory("src/core")
 add_dependencies(core libudunits2)
 add_subdirectory("src/geojson")
+
+if(NGEN_WITH_SQLITE3)
+    add_subdirectory("src/geopackage")
+    target_link_libraries(ngen PUBLIC NGen::geopackage)
+endif()
+
 add_subdirectory("src/realizations/catchment")
 add_subdirectory("src/models/tshirt")
 add_subdirectory("src/models/kernels/reservoir")
diff --git a/include/geopackage/GeoPackage.hpp b/include/geopackage/GeoPackage.hpp
new file mode 100644
index 0000000000..71203f4bcb
--- /dev/null
+++ b/include/geopackage/GeoPackage.hpp
@@ -0,0 +1,62 @@
+#ifndef NGEN_GEOPACKAGE_H
+#define NGEN_GEOPACKAGE_H
+
+#include "FeatureCollection.hpp"
+#include "NGen_SQLite.hpp"
+
+namespace geopackage {
+
+/**
+ * Build a geometry object from GeoPackage WKB.
+ * 
+ * @param[in] row SQLite iterator at the row containing a geometry column
+ * @param[in] geom_col Name of geometry column containing GPKG WKB
+ * @param[out] bounding_box Bounding box of the geometry to output
+ * @return geojson::geometry GPKG WKB converted and projected to a boost geometry model
+ */
+geojson::geometry build_geometry(
+    const sqlite_iter& row,
+    const std::string& geom_col,
+    std::vector<double>& bounding_box
+);
+
+/**
+ * Build properties from GeoPackage table columns.
+ * 
+ * @param[in] row SQLite iterator at the row containing the data columns
+ * @param[in] geom_col Name of geometry column containing GPKG WKB to ignore
+ * @return geojson::PropertyMap PropertyMap of properties from the given row
+ */
+geojson::PropertyMap build_properties(
+    const sqlite_iter& row,
+    const std::string& geom_col
+);
+
+/**
+ * Build a feature from a GPKG table row
+ * 
+ * @param[in] row SQLite iterator at the row to build a feature from
+ * @param[in] geom_col Name of geometry column containing GPKG WKB
+ * @return geojson::Feature Feature containing geometry and properties from the given row
+ */
+geojson::Feature build_feature(
+    const sqlite_iter& row,
+    const std::string& geom_col
+);
+
+/**
+ * Build a feature collection from a GPKG layer
+ *
+ * @param[in] gpkg_path Path to GPKG file
+ * @param[in] layer Layer name within GPKG file to create a collection from
+ * @param[in] ids optional subset of feature IDs to capture (if empty, the entire layer is converted)
+ * @return std::shared_ptr<geojson::FeatureCollection> 
+ */
+std::shared_ptr<geojson::FeatureCollection> read(
+    const std::string& gpkg_path,
+    const std::string& layer,
+    const std::vector<std::string>& ids
+);
+
+} // namespace geopackage
+#endif // NGEN_GEOPACKAGE_H
diff --git a/include/geopackage/NGen_SQLite.hpp b/include/geopackage/NGen_SQLite.hpp
new file mode 100644
index 0000000000..9aaf55d3a7
--- /dev/null
+++ b/include/geopackage/NGen_SQLite.hpp
@@ -0,0 +1,230 @@
+#ifndef NGEN_GEOPACKAGE_SQLITE_H
+#define NGEN_GEOPACKAGE_SQLITE_H
+
+#include <memory>
+#include <stdexcept>
+#include <vector>
+#include <string>
+
+#include <sqlite3.h>
+
+namespace geopackage {
+
+/**
+ * Deleter used to provide smart pointer support for sqlite3 structs.
+ */
+struct sqlite_deleter
+{
+    void operator()(sqlite3* db) { sqlite3_close_v2(db); }
+    void operator()(sqlite3_stmt* stmt) { sqlite3_finalize(stmt); }
+};
+
+/**
+ * Smart pointer (shared) type for sqlite3 prepared statements
+ */
+using stmt_t = std::unique_ptr<sqlite3_stmt, sqlite_deleter>;
+
+/**
+ * SQLite3 row iterator
+ *
+ * Provides a simple iterator-like implementation
+ * over rows of a SQLite3 query.
+ */
+class sqlite_iter
+{
+  private:
+    stmt_t stmt;
+    int iteration_step = -1;
+    bool iteration_finished = false;
+
+    int column_count;
+    std::vector<std::string> column_names;
+
+    // vector of SQLITE data types, see: https://www.sqlite.org/datatype3.html
+    std::vector<int> column_types;
+
+    // returns the raw pointer to the sqlite statement
+    sqlite3_stmt* ptr() const noexcept;
+
+    // checks if int is out of range, and throws error if so
+    void handle_get_index(int) const;
+
+  public:
+    sqlite_iter(stmt_t stmt);
+
+    /**
+     * Check if a row iterator is finished
+     * 
+     * @return true if next() returned SQLITE_DONE
+     * @return false if there is more rows available
+     */
+    bool done() const noexcept;
+
+    /**
+     * Step into the next row of a SQLite query
+     * 
+     * If the query is finished, next() acts idempotently,
+     * but will change done() to return true.
+     * @return sqlite_iter& returns itself
+     */
+    sqlite_iter& next();
+
+    /**
+     * Restart an iteration to its initial state.
+     * next() must be called after calling this.
+     * 
+     * @return sqlite_iter& returns itself
+     */
+    sqlite_iter& restart();
+
+    /**
+     * Get the current row index for the iterator
+     * 
+     * @return int the current row index, or -1 if next() hasn't been called
+     */
+    int current_row() const noexcept;
+
+    /**
+     * Get the number of columns within this iterator
+     * @return int number of columns in query
+     */
+    int num_columns() const noexcept;
+
+    /**
+     * Return the column index for a named column 
+     * 
+     * @param name column name to search for
+     * @return int index of given column name, or -1 if not found.
+     */
+    int column_index(const std::string& name) const noexcept;
+
+    /**
+     * Get a vector of column names
+     * 
+     * @return const std::vector<std::string>& column names as a vector of strings
+     */
+    const std::vector<std::string>& columns() const noexcept { return this->column_names; }
+
+    /**
+     * Get a vector of column types
+     *
+     * See https://www.sqlite.org/datatype3.html for type values. The integers
+     * are the affinity for data types.
+     * @return const std::vector<int>& column types as a vector of ints
+     */
+    const std::vector<int>& types() const noexcept { return this->column_types; }
+
+    /**
+     * Get a column value from a row iterator by index
+     * 
+     * @tparam T Type to parse value as, i.e. int
+     * @param col Column index to parse
+     * @return T value at column `col`
+     */
+    template<typename T>
+    T get(int col) const;
+
+    /**
+     * Get a column value from a row iterator by name
+     * 
+     * @tparam T Type to parse value as, i.e. int
+     * @return T value at the named column
+     */
+    template<typename T>
+    T get(const std::string&) const;
+};
+
+template<typename T>
+inline T sqlite_iter::get(const std::string& name) const
+{
+    const int index = this->column_index(name);
+    return this->get<T>(index);
+}
+
+/**
+ * Wrapper around a SQLite3 database
+ */
+class sqlite
+{
+  private:
+    /**
+     * Smart pointer (unique) type for sqlite3 database
+     */
+    using sqlite_t = std::unique_ptr<sqlite3, sqlite_deleter>;
+  
+    sqlite_t conn = nullptr;
+
+  public:
+    sqlite() = delete;
+
+    /**
+     * Construct a new sqlite object from a path to database
+     *
+     * @param path File path to sqlite3 database
+     */
+    sqlite(const std::string& path);
+
+    sqlite(sqlite& db) = delete;
+  
+    sqlite& operator=(sqlite& db) = delete;
+
+    /**
+     * Take ownership of a sqlite3 database
+     *
+     * @param db sqlite3 database object
+     */
+    sqlite(sqlite&& db) = default;
+
+    /**
+     * Move assignment operator
+     *
+     * @param db sqlite3 database object
+     * @return sqlite& reference to sqlite3 database
+     */
+    sqlite& operator=(sqlite&& db) = default;
+
+    /**
+     * Return the originating sqlite3 database pointer
+     *
+     * @return sqlite3*
+     */
+    sqlite3* connection() const noexcept;
+
+    /**
+     * Check if SQLite database contains a given table
+     *
+     * @param table name of table
+     * @return true if table does exist
+     * @return false if table does not exist
+     */
+    bool has_table(const std::string& table) noexcept;
+
+    /**
+     * Query the SQLite Database and get the result
+     * @param statement String query
+     * @return sqlite_iter SQLite row iterator
+     */
+    sqlite_iter query(const std::string& statement);
+
+    /**
+     * Query the SQLite Database with multiple boundable text parameters.
+     * 
+     * @param statement String query with parameters
+     * @param binds text parameters to bind to statement
+     * @return sqlite_iter SQLite row iterator
+     */
+    sqlite_iter query(const std::string& statement, const std::vector<std::string>& binds);
+
+    /**
+     * Query the SQLite Database with a bound statement and get the result
+     * @param statement String query with parameters
+     * @param params parameters to bind to statement
+     * @return sqlite_iter SQLite row iterator
+     */
+    template<typename... T>
+    sqlite_iter query(const std::string& statement, T const&... params);
+};
+
+} // namespace geopackage
+
+#endif // NGEN_GEOPACKAGE_SQLITE_H
diff --git a/include/geopackage/WKB.hpp b/include/geopackage/WKB.hpp
new file mode 100644
index 0000000000..a4927c017b
--- /dev/null
+++ b/include/geopackage/WKB.hpp
@@ -0,0 +1,149 @@
+#ifndef NGEN_GEOPACKAGE_WKB_H
+#define NGEN_GEOPACKAGE_WKB_H
+
+#include "EndianCopy.hpp"
+#include "JSONGeometry.hpp"
+#include <boost/geometry/srs/epsg.hpp>
+#include <boost/geometry/srs/transformation.hpp>
+
+namespace bg = boost::geometry;
+
+namespace geopackage {
+
+/**
+ * A recursive WKB reader.
+ *
+ * @note This WKB implementation follows a subset of the
+ *       OGC Specification found in https://www.ogc.org/standard/sfa/.
+ *       This is a strict WKB implementation and does not support
+ *       Extended WKB (EWKB) or Tiny WKB (TWKB).
+ */
+struct wkb {
+    using point_t = bg::model::point<double, 2, bg::cs::cartesian>;
+    using linestring_t = bg::model::linestring<point_t>;
+    using polygon_t = bg::model::polygon<point_t>;
+    using multipoint_t = bg::model::multi_point<point_t>;
+    using multilinestring_t = bg::model::multi_linestring<linestring_t>;
+    using multipolygon_t = bg::model::multi_polygon<polygon_t>;
+    using geometry = boost::variant<
+        point_t,
+        linestring_t,
+        polygon_t,
+        multipoint_t,
+        multilinestring_t,
+        multipolygon_t
+    >;
+
+    using byte_t = uint8_t;
+    using byte_vector = std::vector<byte_t>;
+
+    /**
+     * projection visitor. applied with boost to project from
+     * cartesian coordinates to WGS84.
+     */
+    struct wgs84;
+
+    // prevent instatiation of this struct
+    wkb() = delete;
+
+    /**
+     * Read WKB from a given buffer
+     * @param[in] buffer byte vector buffer
+     * @return geometry wkb::geometry struct containing the geometry data from the buffer
+     */
+    static geometry read(const byte_vector& buffer);
+
+    static bg::srs::dpar::parameters<> get_prj(uint32_t srid);
+
+  private:
+    /**
+     * Read a WKB point into a cartesian model.
+     * @return point_t 
+     */
+    static point_t read_point(const byte_vector&, int&, uint8_t);
+
+    /**
+     * Read a WKB linestring into a cartesian model.
+     * @return linestring_t 
+     */
+    static linestring_t read_linestring(const byte_vector&, int&, uint8_t);
+
+    /**
+     * Read a WKB polygon into a cartesian model.
+     * @return polygon_t 
+     */
+    static polygon_t read_polygon(const byte_vector&, int&, uint8_t);
+
+    /**
+     * Read a WKB multipoint into a cartesian model.
+     * @return multipoint_t 
+     */
+    static multipoint_t read_multipoint(const byte_vector&, int&, uint8_t);
+
+    /**
+     * Read a WKB multilinestring into a cartesian model.
+     * @return multilinestring_t 
+     */
+    static multilinestring_t read_multilinestring(const byte_vector&, int&, uint8_t);
+
+    /**
+     * Read a WKB multipolygon into a cartesian model.
+     * @return multipolygon_t 
+     */
+    static multipolygon_t read_multipolygon(const byte_vector&, int&, uint8_t);
+};
+
+/**
+ * EPSG 5070 projection definition for use with boost::geometry.
+ * 
+ * @note this is required because boost 1.72.0 does not
+ *       have an EPSG definition for 5070 in boost::srs::epsg.
+ */
+const auto epsg5070 = bg::srs::dpar::parameters<>(bg::srs::dpar::proj_aea)
+                                                 (bg::srs::dpar::ellps_grs80)
+                                                 (bg::srs::dpar::towgs84, {0,0,0,0,0,0,0})
+                                                 (bg::srs::dpar::lat_0, 23)
+                                                 (bg::srs::dpar::lon_0, -96)
+                                                 (bg::srs::dpar::lat_1, 29.5)
+                                                 (bg::srs::dpar::lat_2, 45.5)
+                                                 (bg::srs::dpar::x_0, 0)
+                                                 (bg::srs::dpar::y_0, 0);
+
+/**
+ * EPSG 3857 projection definition for use with boost::geometry.
+ * 
+ * @note this is required because boost 1.72.0 does not
+ *       have an EPSG definition for 3857 in boost::srs::epsg.
+ */
+const auto epsg3857 = bg::srs::dpar::parameters<>(bg::srs::dpar::proj_merc)
+                                                 (bg::srs::dpar::units_m)
+                                                 (bg::srs::dpar::no_defs)
+                                                 (bg::srs::dpar::a, 6378137)
+                                                 (bg::srs::dpar::b, 6378137)
+                                                 (bg::srs::dpar::lat_ts, 0)
+                                                 (bg::srs::dpar::lon_0, 0)
+                                                 (bg::srs::dpar::x_0, 0)
+                                                 (bg::srs::dpar::y_0, 0)
+                                                 (bg::srs::dpar::k, 1);
+
+struct wkb::wgs84 : public boost::static_visitor<geojson::geometry>
+{
+    wgs84(uint32_t srs, const bg::srs::transformation<>& tr)
+        : srs(srs)
+        , tr(tr) {};
+
+    geojson::geometry operator()(point_t& g);
+    geojson::geometry operator()(linestring_t& g);
+    geojson::geometry operator()(polygon_t& g);
+    geojson::geometry operator()(multipoint_t& g);
+    geojson::geometry operator()(multilinestring_t& g);
+    geojson::geometry operator()(multipolygon_t& g);
+
+    private:
+        uint32_t srs;
+        const bg::srs::transformation<>& tr;
+};
+
+} // namespace geopackage
+
+#endif // NGEN_GEOPACKAGE_WKB_H
diff --git a/include/utilities/EndianCopy.hpp b/include/utilities/EndianCopy.hpp
new file mode 100644
index 0000000000..8bce187052
--- /dev/null
+++ b/include/utilities/EndianCopy.hpp
@@ -0,0 +1,58 @@
+#ifndef NGEN_ENDIANCOPY_H
+#define NGEN_ENDIANCOPY_H
+
+#include <vector>
+
+#include <boost/endian.hpp>
+
+namespace utils {
+
+/**
+ * @brief
+ * Copies bytes from @param{src} to @param{dst},
+ * converts the endianness to native,
+ * and increments @param{index} by the number of bytes
+ * used to store @tparam{S}
+ * 
+ * @tparam S a primitive type
+ * @param src a vector of bytes
+ * @param index an integral type tracking the starting position of @param{dst}'s memory
+ * @param dst output primitive
+ * @param order endianness value (0x01 == Little; 0x00 == Big)
+ */
+template<typename S>
+void copy_from(const std::vector<uint8_t>& src, int& index, S& dst, uint8_t order)
+{
+    std::memcpy(&dst, src.data() + index, sizeof(dst));
+
+    if (order == 0x01) {
+        boost::endian::little_to_native_inplace(dst);
+    } else {
+        boost::endian::big_to_native_inplace(dst);
+    }
+
+    index += sizeof(S);
+}
+
+// boost::endian doesn't support using primitive doubles
+// see: https://github.com/boostorg/endian/issues/36
+template<>
+inline void copy_from<double>(const std::vector<uint8_t>& src, int& index, double& dst, uint8_t order)
+{
+    static_assert(sizeof(uint64_t) == sizeof(double), "sizeof(uint64_t) is not the same as sizeof(double)!");
+
+    uint64_t tmp;
+
+    // copy into uint64_t
+    copy_from(src, index, tmp, order); 
+
+    // copy resolved endianness into double
+    std::memcpy(&dst, &tmp, sizeof(double));
+
+    // above call to copy_from handles index
+    // incrementing, so we don't need to.
+}
+
+}
+
+#endif // NGEN_ENDIANCOPY_H
diff --git a/src/NGen.cpp b/src/NGen.cpp
index 5734820960..07e5b5827c 100644
--- a/src/NGen.cpp
+++ b/src/NGen.cpp
@@ -7,6 +7,10 @@
 #include <Catchment_Formulation.hpp>
 #include <HY_Features.hpp>
 
+#ifdef NGEN_WITH_SQLITE3
+#include <GeoPackage.hpp>
+#endif
+
 #include "NGenConfig.h"
 #include "tshirt_params.h"
 
@@ -259,11 +263,29 @@ int main(int argc, char *argv[]) {
     #endif // NGEN_MPI_ACTIVE
 
     // TODO: Instead of iterating through a collection of FeatureBase objects mapping to nexi, we instead want to iterate through HY_HydroLocation objects
-    geojson::GeoJSON nexus_collection = geojson::read(nexusDataFile, nexus_subset_ids);
+    geojson::GeoJSON nexus_collection;
+    if (boost::algorithm::ends_with(nexusDataFile, "gpkg")) {
+      #ifdef NGEN_WITH_SQLITE3
+      nexus_collection = geopackage::read(nexusDataFile, "nexus", nexus_subset_ids);
+      #else
+      throw std::runtime_error("SQLite3 support required to read GeoPackage files.");
+      #endif
+    } else {
+      nexus_collection = geojson::read(nexusDataFile, nexus_subset_ids);
+    }
     std::cout << "Building Catchment collection" << std::endl;
 
     // TODO: Instead of iterating through a collection of FeatureBase objects mapping to catchments, we instead want to iterate through HY_Catchment objects
-    geojson::GeoJSON catchment_collection = geojson::read(catchmentDataFile, catchment_subset_ids);
+    geojson::GeoJSON catchment_collection;
+    if (boost::algorithm::ends_with(catchmentDataFile, "gpkg")) {
+      #ifdef NGEN_WITH_SQLITE3
+      catchment_collection = geopackage::read(catchmentDataFile, "divides", catchment_subset_ids);
+      #else
+      throw std::runtime_error("SQLite3 support required to read GeoPackage files.");
+      #endif
+    } else {
+      catchment_collection = geojson::read(catchmentDataFile, catchment_subset_ids);
+    }
     
     for(auto& feature: *catchment_collection)
     {
diff --git a/src/geopackage/CMakeLists.txt b/src/geopackage/CMakeLists.txt
new file mode 100644
index 0000000000..acc592bac1
--- /dev/null
+++ b/src/geopackage/CMakeLists.txt
@@ -0,0 +1,26 @@
+cmake_minimum_required(VERSION 3.10)
+
+
+string(COMPARE EQUAL "${CMAKE_CXX_COMPILER_ID}" "IntelLLVM" _cmp)
+if (NOT _cmp)
+    message(WARNING "[NGen::geopackage] wkb.cpp cannot be optimized with "
+                    "${CMAKE_CXX_COMPILER_ID} due to a suspected optimizer/boost issue:\n"
+                    "https://github.com/NOAA-OWP/ngen/issues/567.\n"
+                    "Use IntelLLVM if optimization for this source file is required.")
+    # !! Required due to optimizer issue with either clang or
+    # !! boost::geometry::srs::transformation
+    set_source_files_properties(wkb.cpp PROPERTIES COMPILE_FLAGS -O0)
+endif()
+
+add_library(geopackage geometry.cpp
+                       properties.cpp
+                       feature.cpp
+                       read.cpp
+                       wkb.cpp
+                       sqlite/iterator.cpp
+                       sqlite/database.cpp
+)
+add_library(NGen::geopackage ALIAS geopackage)
+target_include_directories(geopackage PUBLIC ${PROJECT_SOURCE_DIR}/include/geopackage)
+target_include_directories(geopackage PUBLIC ${PROJECT_SOURCE_DIR}/include/utilities)
+target_link_libraries(geopackage PUBLIC NGen::geojson Boost::boost sqlite3)
diff --git a/src/geopackage/feature.cpp b/src/geopackage/feature.cpp
new file mode 100644
index 0000000000..8b51e4d495
--- /dev/null
+++ b/src/geopackage/feature.cpp
@@ -0,0 +1,80 @@
+#include "GeoPackage.hpp"
+
+// Points don't have a bounding box, so we can say its bbox is itself    
+inline void build_point_bbox(const geojson::geometry& geom, std::vector<double>& bbox)
+{
+    const auto& pt = boost::get<geojson::coordinate_t>(geom);
+    bbox[0] = pt.get<0>();
+    bbox[1] = pt.get<1>();
+    bbox[2] = pt.get<0>();
+    bbox[3] = pt.get<1>();
+}
+
+geojson::Feature geopackage::build_feature(
+  const sqlite_iter& row,
+  const std::string& geom_col
+)
+{
+    std::vector<double> bounding_box(4);
+    const auto id                    = row.get<std::string>("id");
+    geojson::PropertyMap properties  = build_properties(row, geom_col);
+    geojson::geometry geometry       = build_geometry(row, geom_col, bounding_box);
+
+    // Convert variant type (0-based) to FeatureType
+    const auto wkb_type = geojson::FeatureType(geometry.which() + 1);
+
+    switch(wkb_type) {
+        case geojson::FeatureType::Point:
+            build_point_bbox(geometry, bounding_box);
+            return std::make_shared<geojson::PointFeature>(
+                boost::get<geojson::coordinate_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::LineString:
+            return std::make_shared<geojson::LineStringFeature>(
+                boost::get<geojson::linestring_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::Polygon:
+            return std::make_shared<geojson::PolygonFeature>(
+                boost::get<geojson::polygon_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::MultiPoint:
+            return std::make_shared<geojson::MultiPointFeature>(
+                boost::get<geojson::multipoint_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::MultiLineString:
+            return std::make_shared<geojson::MultiLineStringFeature>(
+                boost::get<geojson::multilinestring_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::MultiPolygon:
+            return std::make_shared<geojson::MultiPolygonFeature>(
+                boost::get<geojson::multipolygon_t>(geometry),
+                id,
+                properties,
+                bounding_box
+            );
+        case geojson::FeatureType::GeometryCollection:
+            return std::make_shared<geojson::CollectionFeature>(
+                std::vector<geojson::geometry>{geometry},
+                id,
+                properties,
+                bounding_box
+            );
+        default:
+            throw std::runtime_error("invalid WKB feature type. Received: " + std::to_string(geometry.which() + 1));
+    }
+}
diff --git a/src/geopackage/geometry.cpp b/src/geopackage/geometry.cpp
new file mode 100644
index 0000000000..f304b761e8
--- /dev/null
+++ b/src/geopackage/geometry.cpp
@@ -0,0 +1,85 @@
+#include "GeoPackage.hpp"
+#include "EndianCopy.hpp"
+#include "WKB.hpp"
+
+geojson::geometry geopackage::build_geometry(
+    const sqlite_iter& row,
+    const std::string& geom_col,
+    std::vector<double>& bounding_box
+)
+{
+    const std::vector<uint8_t> geometry_blob = row.get<std::vector<uint8_t>>(geom_col);
+    if (geometry_blob[0] != 'G' || geometry_blob[1] != 'P') {
+        throw std::runtime_error("expected geopackage WKB, but found invalid format instead");
+    }
+
+    int index = 3; // skip version
+
+    // flags
+    const bool is_extended  =  geometry_blob[index] & 0x00100000;
+    const bool is_empty     =  geometry_blob[index] & 0x00010000;
+    const uint8_t indicator = (geometry_blob[index] >> 1) & 0x00000111;
+    const uint8_t endian    =  geometry_blob[index] & 0x00000001;
+    index++;
+
+    // Read srs_id
+    uint32_t srs_id;
+    utils::copy_from(geometry_blob, index, srs_id, endian);
+    
+    const auto epsg = wkb::get_prj(srs_id);
+    const bg::srs::transformation<> prj{epsg, wkb::get_prj(4326)};
+    wkb::wgs84 pvisitor{srs_id, prj};
+    
+    if (indicator > 0 & indicator < 5) {
+        // not an empty envelope
+
+        double min_x = 0, max_x = 0, min_y = 0, max_y = 0;
+        utils::copy_from(geometry_blob, index, min_x, endian); // min_x
+        utils::copy_from(geometry_blob, index, max_x, endian); // max_x
+        utils::copy_from(geometry_blob, index, min_y, endian); // min_y
+        utils::copy_from(geometry_blob, index, max_y, endian); // max_y
+
+        // we need to transform the bounding box from its initial SRS
+        // to EPSG: 4326 -- so, we construct a temporary WKB linestring_t type
+        // which will get projected to a geojson::geometry (aka geojson::linestring_t) type.
+
+        // create a wkb::linestring_t bbox object
+        wkb::point_t max{max_x, max_y};
+        wkb::point_t min{min_x, min_y};
+        geojson::coordinate_t max_prj{};
+        geojson::coordinate_t min_prj{};
+
+        // project the raw bounding box
+        if (srs_id == 4326) {
+            max_prj = geojson::coordinate_t{max.get<0>(), max.get<1>()};
+            min_prj = geojson::coordinate_t{min.get<0>(), min.get<1>()};
+        } else {
+            prj.forward(max, max_prj);
+            prj.forward(min, min_prj);
+        }
+
+        // assign the projected values to the bounding_box parameter
+        bounding_box.clear();
+        bounding_box.resize(4); // only 4, not supporting Z or M dims
+        bounding_box[0] = min_prj.get<0>(); // min_x
+        bounding_box[1] = min_prj.get<1>(); // min_y
+        bounding_box[2] = max_prj.get<0>(); // max_x
+        bounding_box[3] = max_prj.get<1>(); // max_y
+
+        // ensure `index` is at beginning of data
+        if (indicator == 2 || indicator == 3) {
+            index += 2 * sizeof(double);
+        } else if (indicator == 4) {
+            index += 4 * sizeof(double);
+        }
+    }
+
+    if (!is_empty) {
+        const std::vector<uint8_t> geometry_data(geometry_blob.begin() + index, geometry_blob.end());
+        auto wkb_geometry = wkb::read(geometry_data);
+        geojson::geometry geometry = boost::apply_visitor(pvisitor, wkb_geometry);
+        return geometry;
+    } else {
+        return geojson::geometry{};
+    }
+}
diff --git a/src/geopackage/properties.cpp b/src/geopackage/properties.cpp
new file mode 100644
index 0000000000..a124ebd55b
--- /dev/null
+++ b/src/geopackage/properties.cpp
@@ -0,0 +1,52 @@
+#include "GeoPackage.hpp"
+#include "JSONProperty.hpp"
+
+geojson::JSONProperty get_property(const geopackage::sqlite_iter& row, const std::string& name, int type)
+{
+    if (type == SQLITE_INTEGER) {
+        auto val = row.get<int>(name);
+        return geojson::JSONProperty(name, val);
+    } else if (type == SQLITE_FLOAT) {
+        auto val = row.get<double>(name);
+        return geojson::JSONProperty(name, val);
+    } else if (type == SQLITE_TEXT) {
+        auto val = row.get<std::string>(name);
+        return geojson::JSONProperty(name, val);
+    } else {
+        return geojson::JSONProperty(name, "null");
+    }
+}
+
+geojson::PropertyMap geopackage::build_properties(
+    const sqlite_iter& row,
+    const std::string& geom_col
+)
+{
+    geojson::PropertyMap properties;
+
+    std::map<std::string, int> property_types;
+    const auto data_cols = row.columns();
+    const auto data_types = row.types();
+
+    std::transform(
+        data_cols.begin(),
+        data_cols.end(),
+        data_types.begin(),
+        std::inserter(property_types, property_types.end()), [](std::string name, int type) {
+            return std::make_pair(name, type);
+        }
+    );
+
+    for (auto& col : property_types) {
+        auto name = col.first;
+        const auto type = col.second;
+        if (name == geom_col) {
+            continue;
+        }
+
+        geojson::JSONProperty property = get_property(row, name, type);
+        properties.emplace(name, std::move(property));
+    }
+
+    return properties;
+}
diff --git a/src/geopackage/read.cpp b/src/geopackage/read.cpp
new file mode 100644
index 0000000000..4cb7dad988
--- /dev/null
+++ b/src/geopackage/read.cpp
@@ -0,0 +1,125 @@
+#include "GeoPackage.hpp"
+
+#include <numeric>
+#include <regex>
+
+void check_table_name(const std::string& table)
+{
+    if (boost::algorithm::starts_with(table, "sqlite_"))
+        throw std::runtime_error("table `" + table + "` is not queryable");
+
+    std::regex allowed("[^-A-Za-z0-9_ ]+");
+    if (std::regex_match(table, allowed))
+        throw std::runtime_error("table `" + table + "` contains invalid characters");
+}
+
+std::shared_ptr<geojson::FeatureCollection> geopackage::read(
+    const std::string& gpkg_path,
+    const std::string& layer = "",
+    const std::vector<std::string>& ids = {}
+)
+{
+    // Check for malicious/invalid layer input
+    check_table_name(layer);
+
+    sqlite db(gpkg_path);
+
+    // Check if layer exists
+    if (!db.has_table(layer)) {
+        // Since the layer doesn't exist, we need to output some additional
+        // debug information with the error. In this case, we add ALL the tables
+        // available in the GPKG, so that if the user sees this error, then it
+        // might've been either a typo or a bad data input, and they can correct.
+        std::string errmsg = "[" + std::string(sqlite3_errmsg(db.connection())) + "] " +
+                             "table " + layer + " does not exist.\n\tTables: ";
+
+        auto errquery = db.query("SELECT name FROM sqlite_master WHERE type='table'");
+        errquery.next();
+        while(!errquery.done()) {
+            errmsg += errquery.get<std::string>(0);
+            errmsg += ", ";
+            errquery.next();
+        }
+
+        throw std::runtime_error(errmsg);
+    }
+
+    // Layer exists, getting statement for it
+    //
+    // this creates a string in the form:
+    //     WHERE id IN (?, ?, ?, ...)
+    // so that it can be bound by SQLite.
+    // This is safer than trying to concatenate
+    // the IDs together.
+    std::string joined_ids = "";
+    if (!ids.empty()) {
+        joined_ids = " WHERE id IN (?";
+        for (size_t i = 1; i < ids.size(); i++) {
+            joined_ids += ", ?";
+        }
+        joined_ids += ")";
+    }
+
+    // Get number of features
+    sqlite_iter query_get_layer_count = db.query("SELECT COUNT(*) FROM " + layer + joined_ids, ids);
+    query_get_layer_count.next();
+    const int layer_feature_count = query_get_layer_count.get<int>(0);
+
+    #ifndef NGEN_QUIET
+    // output debug info on what is read exactly
+    std::cout << "Reading " << layer_feature_count << " features in layer " << layer;
+    if (!ids.empty()) {
+        std::cout << " (id subset:";
+        for (auto& id : ids) {
+            std::cout << " " << id;
+        }
+        std::cout << ")";
+    }
+    std::cout << std::endl;
+    #endif
+
+    // Get layer feature metadata (geometry column name + type)
+    sqlite_iter query_get_layer_geom_meta = db.query("SELECT column_name FROM gpkg_geometry_columns WHERE table_name = ?", layer);
+    query_get_layer_geom_meta.next();
+    const std::string layer_geometry_column = query_get_layer_geom_meta.get<std::string>(0);
+
+    // Get layer
+    sqlite_iter query_get_layer = db.query("SELECT * FROM " + layer + joined_ids, ids);
+    query_get_layer.next();
+
+    // build features out of layer query
+    std::vector<geojson::Feature> features;
+    features.reserve(layer_feature_count);
+    while(!query_get_layer.done()) {
+        geojson::Feature feature = build_feature(
+            query_get_layer,
+            layer_geometry_column
+        );
+
+        features.push_back(feature);
+        query_get_layer.next();
+    }
+
+    // get layer bounding box from features
+    //
+    // GeoPackage contains a bounding box in the SQLite DB,
+    // however, it is in the SRS of the GPKG. By creating
+    // the bbox after the features are built, the projection
+    // is already done. This also should be fairly cheap to do.
+    double min_x = std::numeric_limits<double>::infinity();
+    double min_y = std::numeric_limits<double>::infinity();
+    double max_x = -std::numeric_limits<double>::infinity();
+    double max_y = -std::numeric_limits<double>::infinity();
+    for (const auto& feature : features) {
+        const auto& bbox = feature->get_bounding_box();
+        min_x = bbox[0] < min_x ? bbox[0] : min_x;
+        min_y = bbox[1] < min_y ? bbox[1] : min_y;
+        max_x = bbox[2] > max_x ? bbox[2] : max_x;
+        max_y = bbox[3] > max_y ? bbox[3] : max_y;
+    }
+
+    return std::make_shared<geojson::FeatureCollection>(
+        std::move(features),
+        std::vector<double>({min_x, min_y, max_x, max_y})
+    );
+}
diff --git a/src/geopackage/sqlite/database.cpp b/src/geopackage/sqlite/database.cpp
new file mode 100644
index 0000000000..41a7d418a4
--- /dev/null
+++ b/src/geopackage/sqlite/database.cpp
@@ -0,0 +1,108 @@
+#include "NGen_SQLite.hpp"
+
+using namespace geopackage;
+
+/**
+ * Get a runtime error based on a function and code.
+ *
+ * @param f String denoting the function where the error originated
+ * @param code sqlite3 result code
+ * @param extra additional messages to add to the end of the error
+ * @return std::runtime_error
+ */
+std::runtime_error sqlite_error(const std::string& f, int code, const std::string& extra = "")
+{
+    std::string errmsg = f + " returned code "
+                           + std::to_string(code)
+                           + " (msg: "
+                           + std::string(sqlite3_errstr(code))
+                           + ")";
+
+    if (!extra.empty()) {
+      errmsg += " ";
+      errmsg += extra;
+    }
+
+    return std::runtime_error(errmsg);
+}
+
+sqlite::sqlite(const std::string& path)
+{
+    sqlite3* conn;
+    int      code = sqlite3_open_v2(path.c_str(), &conn, SQLITE_OPEN_READONLY, NULL);
+    if (code != SQLITE_OK) {
+        throw sqlite_error("sqlite3_open_v2", code);
+    }
+    this->conn = sqlite_t(conn);
+}
+
+sqlite3* sqlite::connection() const noexcept
+{
+    return this->conn.get();
+}
+
+bool sqlite::has_table(const std::string& table) noexcept
+{
+    auto q = this->query("SELECT EXISTS(SELECT 1 from sqlite_master WHERE type='table' AND name=?)", table);
+    q.next();
+    return q.get<int>(0);
+};
+
+sqlite_iter sqlite::query(const std::string& statement)
+{
+    sqlite3_stmt* stmt;
+    const auto    cstmt = statement.c_str();
+    const int     code  = sqlite3_prepare_v2(this->connection(), cstmt, statement.length() + 1, &stmt, NULL);
+
+    if (code != SQLITE_OK) {
+        // something happened, can probably switch on result codes
+        // https://www.sqlite.org/rescode.html
+        throw sqlite_error("sqlite3_prepare_v2", code);
+    }
+
+    return sqlite_iter(std::move(stmt_t(stmt)));
+}
+
+sqlite_iter sqlite::query(const std::string& statement, const std::vector<std::string>& binds)
+{
+    sqlite3_stmt* stmt;
+    const auto    cstmt = statement.c_str();
+    const int     code  = sqlite3_prepare_v2(this->connection(), cstmt, statement.length() + 1, &stmt, NULL);
+
+    if (code != SQLITE_OK) {
+        throw sqlite_error("sqlite3_prepare_v2", code);
+    }
+
+    if (!binds.empty()) {
+        for (size_t i = 0; i < binds.size(); i++) {
+            const int code = sqlite3_bind_text(stmt, i + 1, binds[i].c_str(), -1, SQLITE_TRANSIENT);
+            if (code != SQLITE_OK) {
+                throw sqlite_error("sqlite3_bind_text", code);
+            }
+        }
+    }
+
+    return sqlite_iter(std::move(stmt_t(stmt)));
+}
+
+template<typename... T>
+inline sqlite_iter sqlite::query(const std::string& statement, T const&... params)
+{
+    sqlite3_stmt* stmt;
+    const auto    cstmt = statement.c_str();
+    const int     code  = sqlite3_prepare_v2(this->connection(), cstmt, statement.length() + 1, &stmt, NULL);
+
+    if (code != SQLITE_OK) {
+        throw sqlite_error("sqlite3_prepare_v2", code);
+    }
+
+    std::vector<std::string> binds{ { params... } };
+    for (size_t i = 0; i < binds.size(); i++) {
+        const int code = sqlite3_bind_text(stmt, i + 1, binds[i].c_str(), -1, SQLITE_TRANSIENT);
+        if (code != SQLITE_OK) {
+            throw sqlite_error("sqlite3_bind_text", code);
+        }
+    }
+
+    return sqlite_iter(std::move(stmt_t(stmt)));
+}
diff --git a/src/geopackage/sqlite/iterator.cpp b/src/geopackage/sqlite/iterator.cpp
new file mode 100644
index 0000000000..4cd322533f
--- /dev/null
+++ b/src/geopackage/sqlite/iterator.cpp
@@ -0,0 +1,146 @@
+#include <algorithm>
+#include <memory>
+
+#include "NGen_SQLite.hpp"
+
+using namespace geopackage;
+
+// Defined in database.cpp
+extern std::runtime_error sqlite_error(const std::string& f, int code, const std::string& extra = "");
+
+/**
+ * Runtime error for iterations that haven't started
+ */
+const auto sqlite_get_notstarted_error = std::runtime_error(
+    "sqlite iteration is has not started, get() is not callable (call sqlite_iter::next() before)"
+);
+
+/**
+ * Runtime error for iterations that are finished
+ */
+const auto sqlite_get_done_error = std::runtime_error(
+    "sqlite iteration is done, get() is not callable"
+);
+
+sqlite_iter::sqlite_iter(stmt_t stmt)
+    : stmt(std::move(stmt))
+{
+    // sqlite3_column_type requires the last result code to be
+    // SQLITE_ROW, so we need to iterate on the first row.
+    // TODO: need to test how this functions if the last result code
+    //       was SQLITE_DONE.
+    this->next();
+    this->column_count = sqlite3_column_count(this->ptr());
+    this->column_names.reserve(this->column_count);
+    this->column_types.reserve(this->column_count);
+
+    for (int i = 0; i < this->column_count; i++) {
+        this->column_names.emplace_back(sqlite3_column_name(this->ptr(), i));
+        this->column_types.emplace_back(sqlite3_column_type(this->ptr(), i));
+    }
+
+    this->restart();
+}
+
+sqlite3_stmt* sqlite_iter::ptr() const noexcept
+{
+    return this->stmt.get();
+}
+
+void sqlite_iter::handle_get_index(int col) const
+{
+    if (this->done()) {
+        throw sqlite_get_done_error;
+    }
+
+    if (this->current_row() == -1) {
+        throw sqlite_get_notstarted_error;
+    }
+
+    if (col < 0 || col >= this->column_count) {
+        
+        throw std::out_of_range(
+            "column " + std::to_string(col) + " out of range of " + std::to_string(this->column_count) + " columns"
+        );
+    }
+}
+
+bool sqlite_iter::done() const noexcept
+{
+    return this->iteration_finished;
+}
+
+sqlite_iter& sqlite_iter::next()
+{
+    if (!this->done()) {
+        const int returncode = sqlite3_step(this->ptr());
+        if (returncode == SQLITE_DONE) {
+            this->iteration_finished = true;
+        } else if (returncode != SQLITE_ROW) {
+            throw sqlite_error("sqlite3_step", returncode);
+        }
+        this->iteration_step++;
+    }
+
+    return *this;
+}
+
+sqlite_iter& sqlite_iter::restart()
+{
+    sqlite3_reset(this->ptr());
+    this->iteration_step = -1;
+    this->iteration_finished = false;
+    return *this;
+}
+
+int sqlite_iter::current_row() const noexcept
+{
+    return this->iteration_step;
+}
+
+int sqlite_iter::num_columns() const noexcept
+{
+    return this->column_count;
+}
+
+int sqlite_iter::column_index(const std::string& name) const noexcept
+{
+    const size_t pos =
+      std::distance(this->column_names.begin(), std::find(this->column_names.begin(), this->column_names.end(), name));
+
+    return pos >= this->column_names.size() ? -1 : pos;
+}
+
+template<>
+std::vector<uint8_t> sqlite_iter::get<std::vector<uint8_t>>(int col) const
+{
+    this->handle_get_index(col);
+    int size = sqlite3_column_bytes(this->ptr(), col);
+    const uint8_t* ptr = static_cast<const uint8_t*>(sqlite3_column_blob(this->ptr(), col));
+    std::vector<uint8_t> blob(ptr, ptr+size);
+    return blob;
+}
+
+template<>
+double sqlite_iter::get<double>(int col) const
+{
+    this->handle_get_index(col);
+    return sqlite3_column_double(this->ptr(), col);
+}
+
+template<>
+int sqlite_iter::get<int>(int col) const
+{
+    this->handle_get_index(col);
+    return sqlite3_column_int(this->ptr(), col);
+}
+
+template<>
+std::string sqlite_iter::get<std::string>(int col) const
+{
+    this->handle_get_index(col);
+    // TODO: this won't work with non-ASCII text
+    int size = sqlite3_column_bytes(this->ptr(), col);
+    const unsigned char* ptr = sqlite3_column_text(this->ptr(), col);
+    return std::string(ptr, ptr+size);
+}
diff --git a/src/geopackage/wkb.cpp b/src/geopackage/wkb.cpp
new file mode 100644
index 0000000000..a47974799b
--- /dev/null
+++ b/src/geopackage/wkb.cpp
@@ -0,0 +1,339 @@
+#include "WKB.hpp"
+
+using namespace geopackage;
+
+enum wkb_geom_t {
+    geometry            = 0,
+    point               = 1,
+    linestring          = 2,
+    polygon             = 3,
+    multipoint          = 4,
+    multilinestring     = 5,
+    multipolygon        = 6,
+    geometry_collection = 7
+};
+
+void throw_if_not_type(uint32_t given, wkb_geom_t expected)
+{
+    if (given != expected) {
+        throw std::runtime_error(
+            "expected WKB geometry type " +
+            std::to_string(expected) +
+            ", but received " +
+            std::to_string(given)
+        );
+    }
+}
+
+// ----------------------------------------------------------------------------
+// WKB Readers
+// ----------------------------------------------------------------------------
+
+typename wkb::point_t wkb::read_point(const byte_vector& buffer, int& index, uint8_t order)
+{
+    double x, y;
+    utils::copy_from(buffer, index, x, order);
+    utils::copy_from(buffer, index, y, order);
+    return point_t{x, y};
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::linestring_t wkb::read_linestring(const byte_vector& buffer, int& index, uint8_t order)
+{
+    uint32_t count;
+    utils::copy_from(buffer, index, count, order);
+
+    linestring_t linestring;
+    linestring.resize(count);
+    for (auto& child : linestring) {
+        child = read_point(buffer, index, order);
+    }
+
+    return linestring;
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::polygon_t wkb::read_polygon(const byte_vector& buffer, int& index, uint8_t order)
+{
+    uint32_t count;
+    utils::copy_from(buffer, index, count, order);
+
+    polygon_t polygon;
+    
+    if (count > 1) {
+        // polygons only have 1 outer ring,
+        // so any extra vectors are considered to be
+        // inner rings.
+        polygon.inners().resize(count - 1);
+    }
+
+    auto outer = read_linestring(buffer, index, order);
+    polygon.outer().reserve(outer.size());
+    for (auto& p : outer) {
+        polygon.outer().push_back(p);
+    }
+
+    for (uint32_t i = 0; i < count - 1; i++) {
+        auto inner = read_linestring(buffer, index, order);
+        polygon.inners().at(i).reserve(inner.size());
+        for (auto& p : inner) {
+            polygon.inners().at(i).push_back(p);
+        }
+    }
+
+    return polygon;
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::multipoint_t wkb::read_multipoint(const byte_vector& buffer, int& index, uint8_t order)
+{
+    uint32_t count;
+    utils::copy_from(buffer, index, count, order);
+
+    multipoint_t mp;
+    mp.resize(count);
+
+    for (auto& point : mp) {
+        const byte_t new_order = buffer[index];
+        index++;
+
+        uint32_t type;
+        utils::copy_from(buffer, index, type, new_order);
+        throw_if_not_type(type, wkb_geom_t::point);
+
+        point = read_point(buffer, index, new_order);
+    }
+
+    return mp;
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::multilinestring_t wkb::read_multilinestring(const byte_vector& buffer, int& index, uint8_t order)
+{
+    uint32_t count;
+    utils::copy_from(buffer, index, count, order);
+
+    multilinestring_t ml;
+    ml.resize(count);
+    for (auto& line : ml) {
+        const byte_t new_order = buffer[index];
+        index++;
+
+        uint32_t type;
+        utils::copy_from(buffer, index, type, new_order);
+        throw_if_not_type(type, wkb_geom_t::linestring);
+
+        line = read_linestring(buffer, index, new_order);
+    }
+
+    return ml;
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::multipolygon_t wkb::read_multipolygon(const byte_vector& buffer, int& index, uint8_t order)
+{
+    uint32_t count;
+    utils::copy_from(buffer, index, count, order);
+
+    multipolygon_t mpl;
+    mpl.resize(count);
+    for (auto& polygon : mpl) {
+        const byte_t new_order = buffer[index];
+        index++;
+
+        uint32_t type;
+        utils::copy_from(buffer, index, type, new_order);
+        throw_if_not_type(type, wkb_geom_t::polygon);
+
+        polygon = read_polygon(buffer, index, new_order);
+    }
+
+    return mpl;
+}
+
+// ----------------------------------------------------------------------------
+
+typename wkb::geometry wkb::read(const byte_vector& buffer)
+{
+    if (buffer.size() < 5) {
+        throw std::runtime_error("buffer reached end before encountering WKB");
+    }
+
+    int index = 0;
+    const byte_t order = buffer[index];
+    index++;
+
+    uint32_t type;
+    utils::copy_from(buffer, index, type, order);
+
+    switch(type) {
+        case wkb_geom_t::point:
+            return read_point(buffer, index, order);
+        case wkb_geom_t::linestring:
+            return read_linestring(buffer, index, order);
+        case wkb_geom_t::polygon:
+            return read_polygon(buffer, index, order);
+        case wkb_geom_t::multipoint:
+            return read_multipoint(buffer, index, order);
+        case wkb_geom_t::multilinestring:
+            return read_multilinestring(buffer, index, order);
+        case wkb_geom_t::multipolygon:
+            return read_multipolygon(buffer, index, order);
+        default:
+            throw std::runtime_error(
+                "this reader only implements OGC geometry types 1-6, "
+                "but received type " + std::to_string(type)
+            );
+    }
+}
+
+// ----------------------------------------------------------------------------
+// WKB Projection Visitor
+// ----------------------------------------------------------------------------
+
+bg::srs::dpar::parameters<> wkb::get_prj(uint32_t srid) {
+    switch(srid) {
+        case 5070:
+            return epsg5070;
+        case 3857:
+            return epsg3857;
+        default:
+            return bg::projections::detail::epsg_to_parameters(srid);
+    }
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(point_t& g)
+{
+    if (this->srs == 4326) {
+        return geojson::coordinate_t(g.get<0>(), g.get<1>());
+    }
+
+    geojson::coordinate_t h;
+    this->tr.forward(g, h);
+    return h;
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(linestring_t& g)
+{
+    geojson::linestring_t h;
+
+    if (this->srs == 4326) {
+        h.reserve(g.size());
+        for (auto&& gg : g) {
+            h.emplace_back(
+                std::move(gg.get<0>()),
+                std::move(gg.get<1>())
+            );
+        }
+    } else {
+        this->tr.forward(g, h);
+    }
+    return h;
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(polygon_t& g)
+{
+    geojson::polygon_t h;
+
+    if(this->srs == 4326) {
+        h.outer().reserve(g.outer().size());
+        for (auto&& gg : g.outer()) {
+            h.outer().emplace_back(
+                std::move(gg.get<0>()),
+                std::move(gg.get<1>())
+            );
+        }
+
+        h.inners().resize(g.inners().size());
+        auto&& inner_g = g.inners().begin();
+        auto&& inner_h = h.inners().begin();
+        for (; inner_g != g.inners().end(); inner_g++, inner_h++) {
+            inner_h->reserve(inner_g->size());
+            for (auto&& gg : *inner_g) {
+                inner_h->emplace_back(
+                    std::move(gg.get<0>()),
+                    std::move(gg.get<1>())
+                );
+            }
+        }
+    } else {
+        this->tr.forward(g, h);
+    }
+    return h;
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(multipoint_t& g)
+{
+    geojson::multipoint_t h;
+
+    if (this->srs == 4326) {
+        h.reserve(g.size());
+        for (auto&& gg : g) {
+            h.emplace_back(
+                std::move(gg.get<0>()),
+                std::move(gg.get<1>())
+            );
+        }
+    } else {
+        this->tr.forward(g, h);
+    }
+
+    return h;
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(multilinestring_t& g)
+{
+    geojson::multilinestring_t h;
+
+    if (this->srs == 4326) {
+        h.resize(g.size());
+        auto&& line_g = g.begin();
+        auto&& line_h = h.begin();
+        for(; line_g != g.end(); line_g++, line_h++) {
+            *line_h = std::move(
+                boost::get<geojson::linestring_t>(this->operator()(*line_g))
+            );
+        }
+    } else {
+        this->tr.forward(g, h);
+    }
+
+    return h;
+}
+
+// ----------------------------------------------------------------------------
+
+geojson::geometry wkb::wgs84::operator()(multipolygon_t& g)
+{
+    geojson::multipolygon_t h;
+
+    if (this->srs == 4326) {
+        h.resize(g.size());
+        auto&& polygon_g = g.begin();
+        auto&& polygon_h = h.begin();
+        for (; polygon_g != g.end(); polygon_g++, polygon_h++) {
+            *polygon_h = std::move(
+                boost::get<geojson::polygon_t>(this->operator()(*polygon_g))
+            );
+        }
+    } else {
+        this->tr.forward(g, h);
+    }
+
+    return h;
+}
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 722da840d2..a7934beb44 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -1,5 +1,6 @@
 project(test)
 
+include(CTest)
 add_subdirectory(googletest)
 include_directories(${gtest_SOURCE_DIR}/include ${gtest_SOURCE_DIR})
 include_directories(${gmock_SOURCE_DIR}/include ${gmock_SOURCE_DIR})
@@ -84,6 +85,16 @@ add_test(test_geojson
         NGen::geojson
         )
 
+########################## GeoPackage Unit Tests
+if(NGEN_WITH_SQLITE3)
+add_test(test_geopackage
+        3
+        geopackage/WKB_Test.cpp
+        geopackage/SQLite_Test.cpp
+        geopackage/GeoPackage_Test.cpp
+        NGen::geopackage)
+endif()
+
 ########################## Realization Config Unit Tests
 add_test(test_realization_config
         1
diff --git a/test/data/geopackage/example.gpkg b/test/data/geopackage/example.gpkg
new file mode 100644
index 0000000000..09a5e736e3
Binary files /dev/null and b/test/data/geopackage/example.gpkg differ
diff --git a/test/data/geopackage/example_3857.gpkg b/test/data/geopackage/example_3857.gpkg
new file mode 100644
index 0000000000..d4bd67e5e5
Binary files /dev/null and b/test/data/geopackage/example_3857.gpkg differ
diff --git a/test/geopackage/GeoPackage_Test.cpp b/test/geopackage/GeoPackage_Test.cpp
new file mode 100644
index 0000000000..391e3cd150
--- /dev/null
+++ b/test/geopackage/GeoPackage_Test.cpp
@@ -0,0 +1,102 @@
+#include <boost/geometry/io/wkt/write.hpp>
+#include <gtest/gtest.h>
+
+#include "GeoPackage.hpp"
+#include "FileChecker.h"
+
+class GeoPackage_Test : public ::testing::Test
+{
+  protected:
+    void SetUp() override 
+    {
+        this->path = utils::FileChecker::find_first_readable({
+            "test/data/geopackage/example.gpkg",
+            "../test/data/geopackage/example.gpkg",
+            "../../test/data/geopackage/example.gpkg"
+        });
+
+        if (this->path.empty()) {
+            FAIL() << "can't find test/data/geopackage/example.gpkg";
+        }
+
+        this->path2 = utils::FileChecker::find_first_readable({
+            "test/data/geopackage/example_3857.gpkg",
+            "../test/data/geopackage/example_3857.gpkg",
+            "../../test/data/geopackage/example_3857.gpkg"
+        });
+
+        if (this->path2.empty()) {
+            FAIL() << "can't find test/data/geopackage/example_3857.gpkg";
+        }
+    }
+
+    void TearDown() override {};
+
+    std::string path;
+    std::string path2;
+};
+
+TEST_F(GeoPackage_Test, geopackage_read_test)
+{
+    const auto gpkg = geopackage::read(this->path, "test", {});
+    EXPECT_NE(gpkg->find("First"), -1);
+    EXPECT_NE(gpkg->find("Second"), -1);
+    const auto bbox = gpkg->get_bounding_box();
+    EXPECT_EQ(bbox.size(), 4);
+    EXPECT_EQ(bbox[0], 102.0);
+    EXPECT_EQ(bbox[1], 0.0);
+    EXPECT_EQ(bbox[2], 105.0);
+    EXPECT_EQ(bbox[3], 1.0);
+    EXPECT_EQ(2, gpkg->get_size());
+
+    const auto& first = gpkg->get_feature(0);
+    const auto& third = gpkg->get_feature(2);
+    EXPECT_EQ(first->get_id(), "First");
+
+    const auto& point = boost::get<geojson::coordinate_t>(first->geometry());
+    EXPECT_EQ(point.get<0>(), 102.0);
+    EXPECT_EQ(point.get<1>(), 0.5);
+
+    ASSERT_TRUE(third == nullptr);
+}
+
+TEST_F(GeoPackage_Test, geopackage_idsubset_test)
+{
+    const auto gpkg = geopackage::read(this->path, "test", { "First" });
+    EXPECT_NE(gpkg->find("First"), -1);
+    EXPECT_EQ(gpkg->find("Second"), -1);
+
+    const auto& first = gpkg->get_feature(0);
+    EXPECT_EQ(first->get_id(), "First");
+    const auto& point = boost::get<geojson::coordinate_t>(first->geometry());
+    EXPECT_EQ(point.get<0>(), 102.0);
+    EXPECT_EQ(point.get<1>(), 0.5);
+
+    ASSERT_TRUE(gpkg->get_feature(1) == nullptr);
+}
+
+// this test is essentially the same as the above, however, the coordinates
+// are stored in EPSG:3857. When read in, they should convert to EPSG:4326.
+TEST_F(GeoPackage_Test, geopackage_projection_test)
+{
+    const auto gpkg = geopackage::read(this->path2, "example_3857", {});
+    EXPECT_NE(gpkg->find("First"), -1);
+    EXPECT_NE(gpkg->find("Second"), -1);
+    const auto bbox = gpkg->get_bounding_box();
+    EXPECT_EQ(bbox.size(), 4);
+    EXPECT_NEAR(bbox[0], 102.0, 0.0001);
+    EXPECT_NEAR(bbox[1], 0.0, 0.0001);
+    EXPECT_NEAR(bbox[2], 105.0, 0.0001);
+    EXPECT_NEAR(bbox[3], 1.0, 0.0001);
+    EXPECT_EQ(2, gpkg->get_size());
+
+    const auto& first = gpkg->get_feature(0);
+    const auto& third = gpkg->get_feature(2);
+    EXPECT_EQ(first->get_id(), "First");
+
+    const auto& point = boost::get<geojson::coordinate_t>(first->geometry());
+    EXPECT_NEAR(point.get<0>(), 102.0, 0.0001);
+    EXPECT_NEAR(point.get<1>(), 0.5, 0.0001);
+
+    ASSERT_TRUE(third == nullptr);
+}
diff --git a/test/geopackage/SQLite_Test.cpp b/test/geopackage/SQLite_Test.cpp
new file mode 100644
index 0000000000..6bba79c2cb
--- /dev/null
+++ b/test/geopackage/SQLite_Test.cpp
@@ -0,0 +1,111 @@
+#include <gtest/gtest.h>
+
+#include "NGen_SQLite.hpp"
+#include "FileChecker.h"
+
+using namespace geopackage;
+
+class SQLite_Test : public ::testing::Test
+{
+  protected:
+    void SetUp() override
+    {
+        this->path = utils::FileChecker::find_first_readable({
+            "test/data/routing/gauge_01073000.gpkg",
+            "../test/data/routing/gauge_01073000.gpkg",
+            "../../test/data/routing/gauge_01073000.gpkg"
+        });
+
+        if (this->path.empty()) {
+            FAIL() << "can't find gauge_01073000.gpkg";
+        }
+    }
+
+    void TearDown() override {};
+    
+    std::string path;
+    
+};
+
+TEST_F(SQLite_Test, sqlite_access_test)
+{
+    sqlite db {this->path};
+    // user wants metadata
+    EXPECT_TRUE(db.has_table("gpkg_contents"));
+    EXPECT_FALSE(db.has_table("some_fake_table"));
+}
+
+TEST_F(SQLite_Test, sqlite_query_test)
+{
+    sqlite db {this->path};
+
+    if (db.connection() == nullptr) {
+        FAIL() << "database is not loaded";
+    }
+
+    // user provides a query
+    const std::string query = "SELECT * FROM gpkg_contents LIMIT 1";
+    sqlite_iter iter = db.query(query);
+
+    EXPECT_EQ(iter.num_columns(), 10);
+    EXPECT_EQ(iter.columns(), std::vector<std::string>({
+        "table_name",
+        "data_type",
+        "identifier",
+        "description",
+        "last_change",
+        "min_x",
+        "min_y",
+        "max_x",
+        "max_y",
+        "srs_id"
+    }));
+
+    // user iterates over row
+    ASSERT_NO_THROW(iter.next());
+
+    // using column indices
+    EXPECT_EQ(iter.get<std::string>(0), "flowpaths");
+    EXPECT_EQ(iter.get<std::string>(1), "features");
+    EXPECT_EQ(iter.get<std::string>(2), "flowpaths");
+    EXPECT_EQ(iter.get<std::string>(3), "");
+    EXPECT_EQ(iter.get<std::string>(4), "2022-10-25T14:33:51.668Z");
+    EXPECT_EQ(iter.get<double>(5), 1995218.564876059);
+    EXPECT_EQ(iter.get<double>(6), 2502240.321178956);
+    EXPECT_EQ(iter.get<double>(7), 2002525.992495368);
+    EXPECT_EQ(iter.get<double>(8), 2508383.058762011);
+    EXPECT_EQ(iter.get<int>(9), 5070);
+
+    // using column_names
+    EXPECT_EQ(iter.get<std::string>("table_name"), "flowpaths");
+    EXPECT_EQ(iter.get<std::string>("data_type"), "features");
+    EXPECT_EQ(iter.get<std::string>("identifier"), "flowpaths");
+    EXPECT_EQ(iter.get<std::string>("description"), "");
+    EXPECT_EQ(iter.get<std::string>("last_change"), "2022-10-25T14:33:51.668Z");
+    EXPECT_EQ(iter.get<double>("min_x"), 1995218.564876059);
+    EXPECT_EQ(iter.get<double>("min_y"), 2502240.321178956);
+    EXPECT_EQ(iter.get<double>("max_x"), 2002525.992495368);
+    EXPECT_EQ(iter.get<double>("max_y"), 2508383.058762011);
+    EXPECT_EQ(iter.get<int>("srs_id"), 5070);
+
+    // reiteration
+    EXPECT_EQ(iter.current_row(), 0);
+    iter.restart();
+    EXPECT_FALSE(iter.done());
+    EXPECT_EQ(iter.current_row(), -1);
+    ASSERT_ANY_THROW(iter.get<int>(0));
+    iter.next();
+    EXPECT_EQ(iter.current_row(), 0);
+    ASSERT_EQ(iter.get<std::string>(0), "flowpaths");
+
+    // finishing
+    iter.next();
+    EXPECT_TRUE(iter.done());
+    EXPECT_EQ(iter.current_row(), 1);
+    ASSERT_ANY_THROW(iter.get<int>(0));
+
+    // next should be idempotent when iteration is done
+    ASSERT_NO_THROW(iter.next());
+    EXPECT_TRUE(iter.done());
+    EXPECT_EQ(iter.current_row(), 1);
+}
diff --git a/test/geopackage/WKB_Test.cpp b/test/geopackage/WKB_Test.cpp
new file mode 100644
index 0000000000..05a6a9b3cf
--- /dev/null
+++ b/test/geopackage/WKB_Test.cpp
@@ -0,0 +1,262 @@
+#include <boost/geometry/core/cs.hpp>
+#include <gtest/gtest.h>
+
+#include <WKB.hpp>
+
+using namespace geopackage;
+
+class WKB_Test : public ::testing::Test
+{
+  protected:
+    WKB_Test() {}
+
+    ~WKB_Test() override {}
+
+    void SetUp() override {
+        this->wkb_bytes["point_little_endian"] = {
+            0x01,
+            0x01, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10, 0x40
+        };
+
+        this->wkb_bytes["point_big_endian"] = {
+            0x00,
+            0x00, 0x00, 0x00, 0x01,
+            0x40, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x40, 0x10, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00
+        };
+
+        // LINESTRING(30 10, 10 30, 40 40)
+        this->wkb_bytes["linestring"] = {
+            0x01,
+            0x02, 0x00, 0x00, 0x00,
+            0x03, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40
+        };
+        
+        // POLYGON((10.689 -25.092, 34.595 -20.170, 38.814 -35.639, 13.502 -39.155, 10.689 -25.092))
+        this->wkb_bytes["polygon"] = {
+            0x01,
+            0x03, 0x00, 0x00, 0x00,
+            0x01, 0x00, 0x00, 0x00,
+            0x05, 0x00, 0x00, 0x00,
+            0x54, 0xE3, 0xA5, 0x9B, 0xC4, 0x60, 0x25, 0x40,
+            0x64, 0x3B, 0xDF, 0x4F, 0x8D, 0x17, 0x39, 0xC0,
+            0x5C, 0x8F, 0xC2, 0xF5, 0x28, 0x4C, 0x41, 0x40,
+            0xEC, 0x51, 0xB8, 0x1E, 0x85, 0x2B, 0x34, 0xC0,
+            0xD5, 0x78, 0xE9, 0x26, 0x31, 0x68, 0x43, 0x40,
+            0x6F, 0x12, 0x83, 0xC0, 0xCA, 0xD1, 0x41, 0xC0,
+            0x1B, 0x2F, 0xDD, 0x24, 0x06, 0x01, 0x2B, 0x40,
+            0xA4, 0x70, 0x3D, 0x0A, 0xD7, 0x93, 0x43, 0xC0,
+            0x54, 0xE3, 0xA5, 0x9B, 0xC4, 0x60, 0x25, 0x40,
+            0x64, 0x3B, 0xDF, 0x4F, 0x8D, 0x17, 0x39, 0xC0
+        };
+
+        // MULTIPOINT(10 40,40 30,20 20,30 10)
+        this->wkb_bytes["multipoint"] = {
+            0x01,
+            0x04, 0x00, 0x00, 0x00,
+            0x04, 0x00, 0x00, 0x00,
+            0x01,
+            0x01, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x01,
+            0x01, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x01,
+            0x01, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x34, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x34, 0x40,
+            0x01,
+            0x01, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40
+        };
+
+        // MULTILINESTRING((10 10,20 20,10 40),(40 40,30 30,40 20,30 10))
+        this->wkb_bytes["multilinestring"] = {
+            0x01,
+            0x05, 0x00, 0x00, 0x00,
+            0x02, 0x00, 0x00, 0x00,
+            0x01,
+            0x02, 0x00, 0x00, 0x00,
+            0x03, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x34, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x34, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x01,
+            0x02, 0x00, 0x00, 0x00,
+            0x04, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x44, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x34, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x3e, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x24, 0x40
+        };
+
+        // POLYGON ((0 0, 5 0, 5 5, 0 5, 0 0), (2 2, 3 2, 3 3, 2 3, 2 2))
+        this->wkb_bytes["polygon_with_holes"] = {
+            0x01,
+            0x03, 0x00, 0x00, 0x00,
+            0x02, 0x00, 0x00, 0x00,
+            0x05, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x14, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,
+            0x05, 0x00, 0x00, 0x00,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40, 
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x08, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40,
+            0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x40
+        };
+    }
+
+    void TearDown() override {}
+    
+    std::map<std::string, std::vector<uint8_t>> wkb_bytes;
+};
+
+TEST_F(WKB_Test, wkb_point_test) // also tests endianness
+{  
+    const auto geom_big    = boost::get<wkb::point_t>(wkb::read(this->wkb_bytes["point_big_endian"]));
+    const auto geom_little = boost::get<wkb::point_t>(wkb::read(this->wkb_bytes["point_little_endian"]));
+    EXPECT_NEAR(geom_big.get<0>(), geom_little.get<0>(), 1e-6);
+    EXPECT_NEAR(geom_big.get<1>(), geom_little.get<1>(), 1e-6);
+    EXPECT_NEAR(geom_big.get<0>(), 2.0, 1e-6);
+    EXPECT_NEAR(geom_big.get<1>(), 4.0, 1e-6);
+    EXPECT_NEAR(geom_big.get<0>(), 2.0, 1e-6);
+    EXPECT_NEAR(geom_big.get<1>(), 4.0, 1e-6);
+}
+
+TEST_F(WKB_Test, wkb_linestring_test)
+{
+    const wkb::geometry geom = wkb::read(this->wkb_bytes["linestring"]);
+    EXPECT_EQ(geom.which() + 1, 2);
+
+    const wkb::linestring_t& line = boost::get<wkb::linestring_t>(geom);
+    const std::vector<std::pair<double, double>> expected_coordinates = {
+        {30, 10}, {10, 30}, {40, 40}
+    };
+
+    for (int i = 0; i < expected_coordinates.size(); i++) {
+        EXPECT_NEAR(line[i].get<0>(), expected_coordinates[i].first, 1e-4);
+        EXPECT_NEAR(line[i].get<1>(), expected_coordinates[i].second, 1e-4);
+    }
+}
+
+TEST_F(WKB_Test, wkb_polygon_test)
+{    
+    const wkb::geometry geom = wkb::read(this->wkb_bytes["polygon"]);
+    EXPECT_EQ(geom.which() + 1, 3); // +1 since variant.which() is 0-based
+
+    const wkb::polygon_t& poly = boost::get<wkb::polygon_t>(geom);
+    const std::vector<std::pair<double, double>> expected_coordinates = {
+        {10.689, -25.092},
+        {34.595, -20.170},
+        {38.814, -35.639},
+        {13.502, -39.155},
+        {10.689, -25.092}
+    };
+
+    for (int i = 0; i < expected_coordinates.size(); i++) {
+        EXPECT_NEAR(poly.outer().at(i).get<0>(), expected_coordinates[i].first, 1e-4);
+        EXPECT_NEAR(poly.outer().at(i).get<1>(), expected_coordinates[i].second, 1e-4);
+    }
+}
+
+TEST_F(WKB_Test, wkb_polygon_with_holes_test)
+{
+    const wkb::geometry geom = wkb::read(this->wkb_bytes["polygon_with_holes"]);
+    EXPECT_EQ(geom.which() + 1, 3);
+
+    const wkb::polygon_t& poly = boost::get<wkb::polygon_t>(geom);
+    ASSERT_EQ(poly.inners().size(), 1);
+
+    const std::vector<std::pair<double, double>> expected_outer = {
+        {0, 0},
+        {5, 0},
+        {5, 5},
+        {0, 5},
+        {0, 0}
+    };
+
+    for (int i = 0; i < expected_outer.size(); i++) {
+        EXPECT_NEAR(poly.outer()[i].get<0>(), expected_outer[i].first, 1e-4);
+        EXPECT_NEAR(poly.outer()[i].get<1>(), expected_outer[i].second, 1e-4);
+    }
+
+    const std::vector<std::pair<double, double>> expected_inner = {
+        {2, 2},
+        {2, 3},
+        {3, 3},
+        {3, 2},
+        {2, 2}
+    };
+
+    for (int i = 0; i < expected_inner.size(); i++) {
+        EXPECT_NEAR(poly.inners().at(0)[i].get<0>(), expected_inner[i].first, 1e-4);
+        EXPECT_NEAR(poly.inners().at(0)[i].get<1>(), expected_inner[i].second, 1e-4);
+    }
+}
+
+TEST_F(WKB_Test, wkb_multipoint_test)
+{
+    const wkb::geometry geom = wkb::read(this->wkb_bytes["multipoint"]);
+    EXPECT_EQ(geom.which() + 1, 4);
+
+    const wkb::multipoint_t& mp = boost::get<wkb::multipoint_t>(geom);
+    const std::vector<std::pair<double, double>> expected_coordinates = {
+        {10, 40}, {40, 30}, {20, 20}, {30, 10}
+    };
+
+    for (int i = 0; i < expected_coordinates.size(); i++) {
+        EXPECT_NEAR(mp[i].get<0>(), expected_coordinates[i].first, 1e-4);
+        EXPECT_NEAR(mp[i].get<1>(), expected_coordinates[i].second, 1e-4);
+    }
+}
+
+TEST_F(WKB_Test, wkb_multilinestring_test)
+{
+    const wkb::geometry geom = wkb::read(this->wkb_bytes["multilinestring"]);
+    EXPECT_EQ(geom.which() + 1, 5);
+
+    const wkb::multilinestring_t& mp = boost::get<wkb::multilinestring_t>(geom);
+    const std::vector<std::vector<std::pair<double, double>>> expected_coordinates = {
+        { {10, 10}, {20, 20}, {10, 40} },
+        { {40, 40}, {30, 30}, {40, 20}, {30, 10} }
+    };
+
+    for (int i = 0; i < expected_coordinates.size(); i++) {
+        for (int j = 0; j < expected_coordinates[i].size(); j++) {
+            EXPECT_NEAR(mp[i][j].get<0>(), expected_coordinates[i][j].first, 1e-4);
+            EXPECT_NEAR(mp[i][j].get<1>(), expected_coordinates[i][j].second, 1e-4);
+        }
+    }
+}