diff --git a/Project.toml b/Project.toml
index fea5879..bc20e42 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CommonDataModel = "1fbeeb36-5f17-413c-809b-666fb144f157"
 CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb"
+DBInterface = "a10d1c49-ce27-4219-8d33-6db1a4562965"
 FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Format = "1fa38f19-a742-5d3f-a2b9-30dd87b9d5f8"
 GRIBDatasets = "82be9cdb-ee19-4151-bdb3-b400788d9abc"
@@ -27,6 +28,7 @@ PlyIO = "42171d58-473b-503a-8d5f-782019eb09ec"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 ReadVTK = "dc215faf-f008-4882-a9f7-a79a826fadc3"
+SQLite = "0aa819cd-b072-5ff4-a722-6bc24af294d9"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
@@ -41,6 +43,7 @@ CSV = "0.10"
 Colors = "0.12, 0.13"
 CommonDataModel = "0.2, 0.3"
 CoordRefSystems = "0.18"
+DBInterface = "2.5"
 FileIO = "1.16"
 Format = "1.3"
 GRIBDatasets = "0.3, 0.4"
@@ -59,6 +62,7 @@ PlyIO = "1.1"
 PrecompileTools = "1.2"
 PrettyTables = "2.2"
 ReadVTK = "0.2"
+SQLite = "1.6"
 Shapefile = "0.13"
 StaticArrays = "1.6"
 Tables = "1.7"
diff --git a/src/GeoIO.jl b/src/GeoIO.jl
index 81ab4b2..15b21a8 100644
--- a/src/GeoIO.jl
+++ b/src/GeoIO.jl
@@ -56,6 +56,10 @@ import ArchGDAL.GDAL
 # PROJJSON CRS
 import JSON3
 
+# SQLite for GeoPackage
+import SQLite
+import DBInterface
+
 # VTK extensions
 const VTKEXTS = [".vtu", ".vtp", ".vtr", ".vts", ".vti"]
 
@@ -72,7 +76,7 @@ const CDMEXTS = [".grib", ".nc"]
 const FORMATS = [
   (extension=".csv", load="CSV.jl", save="CSV.jl"),
   (extension=".geojson", load="GeoJSON.jl", save="GeoJSON.jl"),
-  (extension=".gpkg", load="ArchGDAL.jl", save="ArchGDAL.jl"),
+  (extension=".gpkg", load="GeoIO.jl", save="GeoIO.jl"),
   (extension=".grib", load="GRIBDatasets.jl", save=""),
   (extension=".gslib", load="GslibIO.jl", save="GslibIO.jl"),
   (extension=".jpeg", load="ImageIO.jl", save="ImageIO.jl"),
@@ -127,6 +131,7 @@ include("extra/csv.jl")
 include("extra/gdal.jl")
 include("extra/geotiff.jl")
 include("extra/gis.jl")
+include("extra/gpkg.jl")
 include("extra/img.jl")
 include("extra/msh.jl")
 include("extra/obj.jl")
diff --git a/src/conversion.jl b/src/conversion.jl
index 8ec44d9..f482b6e 100644
--- a/src/conversion.jl
+++ b/src/conversion.jl
@@ -128,5 +128,8 @@ function togeometry(::GI.MultiPolygonTrait, geom, crs)
   Multi([topolygon(g, CRS) for g in GI.getgeom(geom)])
 end
 
+# Handle geometries that are already Meshes.jl types - return as-is
+geom2meshes(geom::Geometry, crs) = geom
+
 geom2meshes(geom, crs=GI.crs(geom)) = geom2meshes(GI.geomtrait(geom), geom, crs)
 geom2meshes(trait, geom, crs) = togeometry(trait, geom, crs)
diff --git a/src/extra/gis.jl b/src/extra/gis.jl
index 90b7d1f..85c5420 100644
--- a/src/extra/gis.jl
+++ b/src/extra/gis.jl
@@ -50,6 +50,8 @@ function giswrite(fname, geotable; warn, kwargs...)
   elseif endswith(fname, ".parquet")
     CRS = crs(domain(geotable))
     GPQ.write(fname, geotable, (:geometry,), projjson(CRS); kwargs...)
+  elseif endswith(fname, ".gpkg")
+    gpkgwrite(fname, geotable; kwargs...)
   else # fallback to GDAL
     agwrite(fname, geotable; kwargs...)
   end
@@ -63,6 +65,8 @@ function gistable(fname; layer, numtype, kwargs...)
     return GJS.read(fname; numbertype=numtype, kwargs...)
   elseif endswith(fname, ".parquet")
     return GPQ.read(fname; kwargs...)
+  elseif endswith(fname, ".gpkg")
+    return gpkgread(fname; layer, kwargs...)
   else # fallback to GDAL
     data = AG.read(fname; kwargs...)
     return AG.getlayer(data, layer - 1)
diff --git a/src/extra/gpkg.jl b/src/extra/gpkg.jl
new file mode 100644
index 0000000..0499adc
--- /dev/null
+++ b/src/extra/gpkg.jl
@@ -0,0 +1,18 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# Shared WKB type constants
+const WKB_POINT = 0x00000001
+const WKB_LINESTRING = 0x00000002
+const WKB_POLYGON = 0x00000003
+const WKB_MULTIPOINT = 0x00000004
+const WKB_MULTILINESTRING = 0x00000005
+const WKB_MULTIPOLYGON = 0x00000006
+
+const WKB_Z = 0x80000000
+const WKB_M = 0x40000000
+const WKB_ZM = 0xC0000000
+
+include("gpkg/read.jl")
+include("gpkg/write.jl")
\ No newline at end of file
diff --git a/src/extra/gpkg/read.jl b/src/extra/gpkg/read.jl
new file mode 100644
index 0000000..f17b4da
--- /dev/null
+++ b/src/extra/gpkg/read.jl
@@ -0,0 +1,413 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+# Main entry point  
+function gpkgread(fname::String; layer::Int=1, kwargs...)
+  db = SQLite.DB(fname)
+  
+  # 1. Read metadata to extract names of tables and CRS type
+  table_name, geom_column, crs_type = gpkgmetadata(db, layer)
+  
+  # 2. Get all column names from the table
+  columns_result = DBInterface.execute(db, "PRAGMA table_info(\"$table_name\")")
+  all_columns = [row.name for row in Tables.rows(columns_result)]
+  
+  # 3. Build vectors for each column by processing row by row
+  
+  # Get first row to determine column types
+  table_result = DBInterface.execute(db, "SELECT * FROM \"$table_name\" LIMIT 1")
+  first_row_state = iterate(Tables.rows(table_result))
+  
+  if isnothing(first_row_state)
+    error("No data found in table $table_name")
+  end
+  
+  first_row = first_row_state[1]
+  
+  # Initialize column vectors with proper types based on first row
+  column_data = Dict{Symbol, Vector}()
+  for col_name in all_columns
+    col_symbol = Symbol(col_name)
+    sample_value = Tables.getcolumn(first_row, col_symbol)
+    if col_name == geom_column
+      column_data[col_symbol] = Any[]  # Geometry column stays as Any for mixed types
+    else
+      # Use the actual type from SQLite for attribute columns
+      column_data[col_symbol] = Vector{typeof(sample_value)}()
+    end
+  end
+  
+  # Process first row
+  for col_name in all_columns
+    col_symbol = Symbol(col_name)
+    if col_name == geom_column
+      # Parse geometry blob
+      blob = Tables.getcolumn(first_row, col_symbol)
+      if ismissing(blob) || isnothing(blob)
+        push!(column_data[col_symbol], blob)
+      else
+        push!(column_data[col_symbol], parse_wkb_geometry_from_blob(blob, crs_type))
+      end
+    else
+      # Store attribute data as-is
+      push!(column_data[col_symbol], Tables.getcolumn(first_row, col_symbol))
+    end
+  end
+  
+  # Query and process remaining rows
+  table_result = DBInterface.execute(db, "SELECT * FROM \"$table_name\"")
+  row_count = 0
+  for row in Tables.rows(table_result)
+    row_count += 1
+    if row_count == 1
+      continue  # Skip first row since we already processed it
+    end
+    
+    for col_name in all_columns
+      col_symbol = Symbol(col_name)
+      if col_name == geom_column
+        # Parse geometry blob
+        blob = Tables.getcolumn(row, col_symbol)
+        if ismissing(blob) || isnothing(blob)
+          push!(column_data[col_symbol], blob)
+        else
+          push!(column_data[col_symbol], parse_wkb_geometry_from_blob(blob, crs_type))
+        end
+      else
+        # Store attribute data as-is
+        push!(column_data[col_symbol], Tables.getcolumn(row, col_symbol))
+      end
+    end
+  end
+  
+  # 4. Build final table
+  final_names = Symbol[]
+  final_columns = Any[]
+  
+  for col_name in all_columns
+    if col_name == geom_column
+      push!(final_names, :geometry)
+    else
+      push!(final_names, Symbol(col_name))
+    end
+    push!(final_columns, column_data[Symbol(col_name)])
+  end
+  
+  close(db)
+  
+  return (; zip(final_names, final_columns)...)
+end
+
+function gpkgmetadata(db::SQLite.DB, layer::Int)
+  # Get feature table names
+  contents_result = DBInterface.execute(db,
+    """
+    SELECT table_name, data_type
+    FROM gpkg_contents
+    WHERE data_type = 'features'
+    ORDER BY table_name
+    """)
+  
+  # Use Iterators.drop to get the requested layer
+  rows = Tables.rows(contents_result)
+  layer_iter = Iterators.drop(rows, layer - 1)
+  
+  if isempty(layer_iter)
+    error("Layer $layer not found in GeoPackage")
+  end
+  first_row = first(layer_iter)
+  
+  table_name = first_row.table_name
+  
+  # Get geometry column information
+  geom_result = DBInterface.execute(db,
+    """
+    SELECT column_name, srs_id
+    FROM gpkg_geometry_columns
+    WHERE table_name = ?
+    """, [table_name])
+  geom_info = first(Tables.rows(geom_result))
+  geom_column = geom_info.column_name
+  srs_id = geom_info.srs_id
+  
+  # Get CRS from SRID according to GeoPackage specification
+  # Handle special cases per GeoPackage spec
+  if srs_id == 0
+    crs_type = LatLon{WGS84Latest}
+  elseif srs_id == -1
+    # Undefined Cartesian coordinate reference system - determine dimensionality from geometry_columns
+    geom_cols_result = DBInterface.execute(db,
+      """
+      SELECT z
+      FROM gpkg_geometry_columns
+      WHERE table_name = ? AND column_name = ?
+      """, [table_name, geom_column])
+    
+    geom_state = iterate(Tables.rows(geom_cols_result))
+    if !isnothing(geom_state)
+      z_flag = geom_state[1].z
+      if z_flag == 1
+        crs_type = Cartesian{NoDatum,3}
+      else
+        crs_type = Cartesian{NoDatum,2}
+      end
+    else
+      crs_type = Cartesian{NoDatum,2}
+    end
+  else
+    # Query the database for the organization and organization_coordsys_id
+    srs_result = DBInterface.execute(db,
+      """
+      SELECT organization, organization_coordsys_id
+      FROM gpkg_spatial_ref_sys
+      WHERE srs_id = ?
+      """, [srs_id])
+    
+    org_state = iterate(Tables.rows(srs_result))
+    if isnothing(org_state)
+      # SRID not found - use fallback to Cartesian (allowed by spec)
+      crs_type = Cartesian{NoDatum}
+    else
+      org_info = org_state[1]
+      
+      # Per GeoPackage spec, organization and organization_coordsys_id are NOT NULL
+      org_name = uppercase(org_info.organization)
+      coord_id = org_info.organization_coordsys_id
+      
+      if org_name == "EPSG"
+        crs_type = CoordRefSystems.get(EPSG{coord_id})
+      elseif org_name == "ESRI"
+        crs_type = CoordRefSystems.get(ESRI{coord_id})
+      else
+        # Per spec: other organizations or NONE
+        crs_type = Cartesian{NoDatum}
+      end
+    end
+  end
+  
+  return table_name, geom_column, crs_type
+end
+
+
+# GeoPackage Binary Header structure
+struct GeoPackageBinaryHeader
+  magic::UInt16
+  version::UInt8
+  flags::UInt8
+  srs_id::Int32
+  min_x::Float64
+  max_x::Float64
+  min_y::Float64
+  max_y::Float64
+  min_z::Union{Float64,Nothing}
+  max_z::Union{Float64,Nothing}
+  min_m::Union{Float64,Nothing}
+  max_m::Union{Float64,Nothing}
+end
+
+function parse_gpb_header(blob::Vector{UInt8})
+  io = IOBuffer(blob)
+  
+  magic = read(io, UInt16)
+  magic == 0x5047 || error("Invalid GeoPackage binary header")
+  
+  version = read(io, UInt8)
+  flags = read(io, UInt8)
+  
+  srs_id = ltoh(read(io, Int32))
+  
+  envelope_type = flags >> 1 & 0x07
+  has_envelope = envelope_type > 0
+  
+  if has_envelope && envelope_type >= 1
+    min_x = ltoh(read(io, Float64))
+    max_x = ltoh(read(io, Float64))
+    min_y = ltoh(read(io, Float64))
+    max_y = ltoh(read(io, Float64))
+  else
+    # Default envelope when not present
+    min_x = max_x = min_y = max_y = 0.0
+  end
+  
+  # Initialize z and m coordinates based on envelope type
+  if has_envelope && envelope_type >= 2
+    min_z = ltoh(read(io, Float64))
+    max_z = ltoh(read(io, Float64))
+  else
+    min_z = max_z = nothing
+  end
+  
+  if has_envelope && envelope_type >= 3
+    min_m = ltoh(read(io, Float64))
+    max_m = ltoh(read(io, Float64))
+  else
+    min_m = max_m = nothing
+  end
+  
+  header = GeoPackageBinaryHeader(
+    magic, version, flags, srs_id,
+    min_x, max_x, min_y, max_y,
+    min_z, max_z, min_m, max_m
+  )
+  
+  return header, position(io)
+end
+
+# Parse GeoPackage Binary (GPB) format - combines header and WKB geometry
+function parse_wkb_geometry_from_blob(blob::Vector{UInt8}, crs_type)
+  header, offset = parse_gpb_header(blob)
+  io = IOBuffer(blob[offset+1:end])
+  parse_wkb_geometry(io, crs_type)
+end
+
+function parse_wkb_geometry(io::IOBuffer, crs_type)
+  byte_order = read(io, UInt8)
+  read_value = byte_order == 0x01 ? ltoh : ntoh
+  
+  wkb_type = read_value(read(io, UInt32))
+  base_type = wkb_type & 0x0FFFFFFF
+  
+  # Determine coordinate dimensionality from CRS type
+  ndims = CoordRefSystems.ncoords(crs_type)
+  
+  # Create parse_coords closure that reads coordinates based on CRS dimensionality
+  parse_coords = io -> read_coords_from_buffer(io, read_value, ndims)
+
+  if base_type == WKB_POINT
+    return parse_wkb_point(io, parse_coords, crs_type)
+  elseif base_type == WKB_LINESTRING
+    return parse_wkb_linestring(io, parse_coords, crs_type)
+  elseif base_type == WKB_POLYGON
+    return parse_wkb_polygon(io, parse_coords, crs_type)
+  elseif base_type == WKB_MULTIPOINT
+    return parse_wkb_multipoint(io, parse_coords, crs_type)
+  elseif base_type == WKB_MULTILINESTRING
+    return parse_wkb_multilinestring(io, parse_coords, crs_type)
+  elseif base_type == WKB_MULTIPOLYGON
+    return parse_wkb_multipolygon(io, parse_coords, crs_type)
+  else
+    error("Unsupported WKB geometry type: $base_type")
+  end
+end
+
+# Read coordinates from buffer based on CRS dimensionality
+function read_coords_from_buffer(io::IOBuffer, read_value::Function, ndims::Int)
+  
+  if ndims == 2
+    x = read_value(read(io, Float64))
+    y = read_value(read(io, Float64))
+    return (x, y)
+  elseif ndims == 3
+    x = read_value(read(io, Float64))
+    y = read_value(read(io, Float64))
+    z = read_value(read(io, Float64))
+    return (x, y, z)
+  else
+    error("Unsupported number of dimensions: $ndims")
+  end
+end
+
+
+# Convert coordinate buffer to CRS coordinates
+function create_crs_coord(buff::Tuple, crs_type::Type{<:CRS})
+  if crs_type <: LatLon
+    crs_type(buff[2], buff[1])  # LatLon(lat, lon) from WKB(x=lon, y=lat)
+  elseif crs_type <: LatLonAlt
+    crs_type(buff[2], buff[1], buff[3])  # LatLonAlt(lat, lon, alt)
+  else
+    crs_type(buff...)  # fallback for other coordinate systems
+  end
+end
+
+function parse_wkb_point(io::IOBuffer, parse_coords::Function, crs_type)
+  buff = parse_coords(io)
+  coord = create_crs_coord(buff, crs_type)
+  return Point(coord)
+end
+
+# Batch parse multiple points efficiently
+function parse_wkb_points(io::IOBuffer, n::Int, parse_coords::Function, crs_type::Type{<:CRS})
+  # Read n coordinate buffers and convert to points
+  map(1:n) do _
+    buff = parse_coords(io)
+    coord = create_crs_coord(buff, crs_type)
+    Point(coord)
+  end
+end
+
+function parse_wkb_linestring(io::IOBuffer, parse_coords::Function, crs_type; is_ring::Bool=false)
+  num_points = Int(ltoh(read(io, UInt32)))
+  
+  # Use batch processing
+  points = parse_wkb_points(io, num_points, parse_coords, crs_type)
+  
+  # Check if WKB writer included duplicate closing point (varies by implementation)
+  has_equal_ends = length(points) >= 2 && points[1] == points[end]
+  
+  # Two-flag logic explained:
+  # - is_ring: Semantic flag indicating this linestring is part of a polygon structure
+  # - has_equal_ends: WKB format flag indicating if writer included duplicate closing point
+  if is_ring || has_equal_ends
+    # Return a Ring, exclude duplicate last point if present
+    if has_equal_ends
+      return Ring(points[1:end-1]...)
+    else
+      return Ring(points...)
+    end
+  else
+    # It's an open linestring, return a Rope
+    return Rope(points...)
+  end
+end
+
+function parse_wkb_polygon(io::IOBuffer, parse_coords::Function, crs_type)
+  num_rings = Int(ltoh(read(io, UInt32)))
+  
+  # Parse exterior ring using ring parsing function
+  exterior_ring = parse_wkb_linestring(io, parse_coords, crs_type; is_ring=true)
+  
+  # Parse interior rings (holes) using map-do syntax
+  hole_rings = map(2:num_rings) do ring_idx
+    parse_wkb_linestring(io, parse_coords, crs_type; is_ring=true)
+  end
+  
+  return PolyArea(exterior_ring, hole_rings...)
+end
+
+
+function parse_wkb_multipoint(io::IOBuffer, parse_coords::Function, crs_type)
+  num_points = Int(ltoh(read(io, UInt32)))
+  
+  points = map(1:num_points) do _
+    # Skip sub-geometry header (byte_order and wkb_type)
+    skip(io, 5)
+    parse_wkb_point(io, parse_coords, crs_type)
+  end
+  
+  return Multi(points)
+end
+
+function parse_wkb_multilinestring(io::IOBuffer, parse_coords::Function, crs_type)
+  num_lines = Int(ltoh(read(io, UInt32)))
+  
+  lines = map(1:num_lines) do _
+    # Skip sub-geometry header (byte_order and wkb_type)
+    skip(io, 5)
+    parse_wkb_linestring(io, parse_coords, crs_type)
+  end
+  
+  return Multi(lines)
+end
+
+function parse_wkb_multipolygon(io::IOBuffer, parse_coords::Function, crs_type)
+  num_polygons = Int(ltoh(read(io, UInt32)))
+  
+  polygons = map(1:num_polygons) do _
+    # Skip sub-geometry header (byte_order and wkb_type)
+    skip(io, 5)
+    parse_wkb_polygon(io, parse_coords, crs_type)
+  end
+  
+  return Multi(polygons)
+end
\ No newline at end of file
diff --git a/src/extra/gpkg/write.jl b/src/extra/gpkg/write.jl
new file mode 100644
index 0000000..7dd5c42
--- /dev/null
+++ b/src/extra/gpkg/write.jl
@@ -0,0 +1,441 @@
+# ------------------------------------------------------------------
+# Licensed under the MIT License. See LICENSE in the project root.
+# ------------------------------------------------------------------
+
+
+function create_gpb_header(srs_id::Int32, geom::Geometry, envelope_type::Int=1)
+  io = IOBuffer()
+  
+  write(io, UInt8(0x47), UInt8(0x50))
+  write(io, UInt8(0))
+  
+  flags = UInt8(0x20 | (envelope_type << 1))
+  write(io, flags)
+  
+  write(io, htol(srs_id))
+  
+  if envelope_type >= 1
+    bbox = boundingbox(geom)
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.min))[1])))
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.max))[1])))
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.min))[2])))
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.max))[2])))
+  end
+  
+  if envelope_type >= 2 && length(CoordRefSystems.raw(coords(bbox.min))) >= 3
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.min))[3])))
+    write(io, htol(Float64(CoordRefSystems.raw(coords(bbox.max))[3])))
+  end
+  
+  return take!(io)
+end
+
+function write_coords(io::IOBuffer, point::Point, dim::Int)
+  coords_pt = CoordRefSystems.raw(coords(point))
+  write(io, Float64(coords_pt[1]))
+  write(io, Float64(coords_pt[2]))
+  if dim >= 3
+    write(io, Float64(coords_pt[3]))
+  end
+end
+
+function write_wkb_point(io::IOBuffer, point::Point)
+  dim = CoordRefSystems.ncoords(typeof(coords(point)))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_POINT
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  write(io, UInt32(wkb_type))
+  
+  write_coords(io, point, dim)
+end
+
+function write_wkb_linestring(io::IOBuffer, rope::Rope)
+  points = vertices(rope)
+  dim = CoordRefSystems.ncoords(typeof(coords(first(points))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_LINESTRING
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  write(io, UInt32(wkb_type))
+  
+  write(io, UInt32(length(points)))
+  
+  for point in points
+    write_coords(io, point, dim)
+  end
+end
+
+function write_wkb_ring_as_linestring(io::IOBuffer, ring::Ring)
+  points = vertices(ring)
+  dim = CoordRefSystems.ncoords(typeof(coords(first(points))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_LINESTRING
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  write(io, UInt32(wkb_type))
+  
+  write(io, UInt32(length(points) + 1))
+  
+  for point in points
+    write_coords(io, point, dim)
+  end
+  
+  write_coords(io, first(points), dim)
+end
+
+function write_wkb_polygon(io::IOBuffer, poly::PolyArea)
+  exterior_ring = boundary(poly)
+  all_rings = [exterior_ring]
+  
+  dim = CoordRefSystems.ncoords(typeof(coords(first(vertices(exterior_ring)))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_POLYGON
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  write(io, UInt32(wkb_type))
+  
+  write(io, UInt32(length(all_rings)))
+  
+  for ring in all_rings
+    points = vertices(ring)
+    write(io, UInt32(length(points) + 1))
+    
+    for point in points
+      write_coords(io, point, dim)
+    end
+    
+    write_coords(io, first(points), dim)
+  end
+end
+
+function write_wkb_multipoint(io::IOBuffer, multi::Multi)
+  parts = parent(multi)
+  dim = CoordRefSystems.ncoords(typeof(coords(first(parts))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_MULTIPOINT
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  
+  write(io, UInt32(wkb_type))
+  write(io, UInt32(length(parts)))
+  
+  for part in parts
+    write_wkb_point(io, part)
+  end
+end
+
+function write_wkb_multilinestring(io::IOBuffer, multi::Multi)
+  parts = parent(multi)
+  dim = CoordRefSystems.ncoords(typeof(coords(first(parts))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_MULTILINESTRING
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  
+  write(io, UInt32(wkb_type))
+  write(io, UInt32(length(parts)))
+  
+  for part in parts
+    write_wkb_linestring(io, part)
+  end
+end
+
+function write_wkb_multipolygon(io::IOBuffer, multi::Multi)
+  parts = parent(multi)
+  dim = CoordRefSystems.ncoords(typeof(coords(first(parts))))
+  
+  write(io, UInt8(0x01))
+  
+  wkb_type = WKB_MULTIPOLYGON
+  if dim == 3
+    wkb_type |= WKB_Z
+  end
+  
+  write(io, UInt32(wkb_type))
+  write(io, UInt32(length(parts)))
+  
+  for part in parts
+    write_wkb_polygon(io, part)
+  end
+end
+
+function write_wkb_geometry(io::IOBuffer, geom::Geometry)
+  if geom isa Point
+    write_wkb_point(io, geom)
+  elseif geom isa Rope
+    write_wkb_linestring(io, geom)
+  elseif geom isa PolyArea
+    write_wkb_polygon(io, geom)
+  elseif geom isa Ring
+    write_wkb_ring_as_linestring(io, geom)
+  elseif geom isa Multi
+    parts = parent(geom)
+    first_part = first(parts)
+    
+    if first_part isa Point
+      write_wkb_multipoint(io, geom)
+    elseif first_part isa Rope
+      write_wkb_multilinestring(io, geom)
+    elseif first_part isa PolyArea
+      write_wkb_multipolygon(io, geom)
+    else
+      error("Unsupported multi-geometry type")
+    end
+  else
+    error("Unsupported geometry type: $(typeof(geom))")
+  end
+end
+
+function create_gpb(geom::Geometry, srs_id::Int32)
+  header = create_gpb_header(srs_id, geom)
+  
+  wkb_io = IOBuffer()
+  write_wkb_geometry(wkb_io, geom)
+  wkb_data = take!(wkb_io)
+  
+  return vcat(header, wkb_data)
+end
+
+function get_srid_for_crs(db::SQLite.DB, crs)
+  if isnothing(crs)
+    return Int32(0)
+  elseif crs <: Cartesian{NoDatum}
+    return Int32(-1)
+  end
+  
+  # Try to extract EPSG code from the CRS type
+  crs_string = string(crs)
+  
+  # Check if it's an EPSG-based CRS
+  if occursin("EPSG", crs_string)
+    # Try to extract the EPSG code
+    m = match(r"EPSG\{(\d+)\}", crs_string)
+    if !isnothing(m)
+      return Int32(parse(Int, m.captures[1]))
+    end
+  end
+  
+  # Handle specific well-known CRS types
+  if crs <: LatLon{WGS84Latest} || crs <: LatLon{WGS84{1762}}
+    return Int32(4326)
+  elseif crs <: Cartesian{WGS84Latest} || crs <: Cartesian{WGS84{1762}}
+    # Use SRID -1 for undefined Cartesian systems
+    return Int32(-1)
+  elseif crs <: Cartesian
+    # Any other Cartesian system
+    return Int32(-1)
+  end
+  
+  return Int32(4326)
+end
+
+function ensure_gpkg_tables(db::SQLite.DB)
+  DBInterface.execute(db, """
+    CREATE TABLE IF NOT EXISTS gpkg_spatial_ref_sys (
+      srs_name TEXT NOT NULL,
+      srs_id INTEGER NOT NULL PRIMARY KEY,
+      organization TEXT NOT NULL,
+      organization_coordsys_id INTEGER NOT NULL,
+      definition TEXT NOT NULL,
+      description TEXT
+    )
+  """)
+  
+  DBInterface.execute(db, """
+    INSERT OR IGNORE INTO gpkg_spatial_ref_sys 
+    (srs_name, srs_id, organization, organization_coordsys_id, definition, description)
+    VALUES 
+    ('WGS 84', 4326, 'EPSG', 4326, 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]', 'WGS 84'),
+    ('Undefined geographic SRS', 0, 'NONE', 0, 'undefined', 'undefined geographic coordinate reference system'),
+    ('Undefined Cartesian SRS', -1, 'NONE', -1, 'undefined', 'undefined Cartesian coordinate reference system')
+  """)
+  
+  DBInterface.execute(db, """
+    CREATE TABLE IF NOT EXISTS gpkg_contents (
+      table_name TEXT NOT NULL PRIMARY KEY,
+      data_type TEXT NOT NULL,
+      identifier TEXT UNIQUE,
+      description TEXT DEFAULT '',
+      last_change DATETIME NOT NULL DEFAULT (strftime('%Y-%m-%dT%H:%M:%fZ','now')),
+      min_x DOUBLE,
+      min_y DOUBLE,
+      max_x DOUBLE,
+      max_y DOUBLE,
+      srs_id INTEGER,
+      CONSTRAINT fk_gc_srs FOREIGN KEY (srs_id) REFERENCES gpkg_spatial_ref_sys(srs_id)
+    )
+  """)
+  
+  DBInterface.execute(db, """
+    CREATE TABLE IF NOT EXISTS gpkg_geometry_columns (
+      table_name TEXT NOT NULL,
+      column_name TEXT NOT NULL,
+      geometry_type_name TEXT NOT NULL,
+      srs_id INTEGER NOT NULL,
+      z TINYINT NOT NULL,
+      m TINYINT NOT NULL,
+      CONSTRAINT pk_geom_cols PRIMARY KEY (table_name, column_name),
+      CONSTRAINT fk_gc_tn FOREIGN KEY (table_name) REFERENCES gpkg_contents(table_name),
+      CONSTRAINT fk_gc_srs FOREIGN KEY (srs_id) REFERENCES gpkg_spatial_ref_sys(srs_id)
+    )
+  """)
+end
+
+function infer_geometry_type(geoms::AbstractVector{<:Geometry})
+  if isempty(geoms)
+    return "GEOMETRY"
+  end
+  
+  # Check for homogeneous geometry types using eltype dispatch
+  T = eltype(geoms)
+  
+  if T <: Point
+    return "POINT"
+  elseif T <: Rope
+    return "LINESTRING"
+  elseif T <: Ring
+    return "LINESTRING"  # Ring is saved as LINESTRING to match GDAL behavior
+  elseif T <: PolyArea
+    return "POLYGON"
+  elseif T <: Multi
+    # Check the element type of Multi
+    element_type = eltype(parent(first(geoms)))
+    if element_type <: Point
+      return "MULTIPOINT"
+    elseif element_type <: Rope
+      return "MULTILINESTRING"
+    elseif element_type <: PolyArea
+      return "MULTIPOLYGON"
+    end
+  end
+  
+  # Fallback for mixed or unknown types
+  return "GEOMETRY"
+end
+
+function gpkgwrite(fname::String, geotable; layer::String="features", kwargs...)
+  db = SQLite.DB(fname)
+  
+  ensure_gpkg_tables(db)
+  
+  # Drop existing table if it exists to avoid data duplication  
+  DBInterface.execute(db, "DROP TABLE IF EXISTS \"$layer\"")
+  # Clean up existing metadata (ignore if tables don't exist)
+  DBInterface.execute(db, "DELETE FROM gpkg_contents WHERE table_name = ? AND EXISTS (SELECT 1 FROM gpkg_contents WHERE table_name = ?)", [layer, layer])
+  DBInterface.execute(db, "DELETE FROM gpkg_geometry_columns WHERE table_name = ? AND EXISTS (SELECT 1 FROM gpkg_geometry_columns WHERE table_name = ?)", [layer, layer])
+  
+  table = values(geotable)
+  domain = GeoTables.domain(geotable)
+  
+  geoms = collect(domain)
+  crs = GeoTables.crs(domain)
+  srs_id = get_srid_for_crs(db, crs)
+  
+  # Handle case where table is Nothing (no attribute columns)
+  cols = isnothing(table) ? Symbol[] : Tables.columnnames(Tables.columns(table))
+  
+  # Start with basic table structure (no automatic fid column)
+  create_table_sql = """
+    CREATE TABLE IF NOT EXISTS "$layer" (
+      geom BLOB NOT NULL
+    """
+  
+  # Add columns only if table is not Nothing
+  if !isnothing(table)
+    for col in cols
+      col_str = string(col)
+      coltype = eltype(Tables.getcolumn(Tables.columns(table), col))
+      
+      if coltype <: Integer
+        sql_type = "INTEGER"
+      elseif coltype <: AbstractFloat
+        sql_type = "REAL"
+      else
+        sql_type = "TEXT"
+      end
+      
+      create_table_sql *= ",\n    \"$col_str\" $sql_type"
+    end
+  end
+  
+  create_table_sql *= "\n)"
+  
+  DBInterface.execute(db, create_table_sql)
+  
+  bbox = boundingbox(domain)
+  min_coords = CoordRefSystems.raw(coords(bbox.min))
+  max_coords = CoordRefSystems.raw(coords(bbox.max))
+  
+  DBInterface.execute(db, """
+    INSERT OR REPLACE INTO gpkg_contents 
+    (table_name, data_type, identifier, min_x, min_y, max_x, max_y, srs_id)
+    VALUES (?, 'features', ?, ?, ?, ?, ?, ?)
+  """, [layer, layer, min_coords[1], min_coords[2], max_coords[1], max_coords[2], srs_id])
+  
+  geom_type = infer_geometry_type(geoms)
+  first_geom = first(geoms)
+  z_flag = if first_geom isa Point
+    CoordRefSystems.ncoords(typeof(coords(first_geom))) >= 3 ? 1 : 0
+  else
+    # For non-Point geometries, get coordinate dimensionality from vertices
+    first_point = first(vertices(first_geom))
+    CoordRefSystems.ncoords(typeof(coords(first_point))) >= 3 ? 1 : 0
+  end
+  
+  DBInterface.execute(db, """
+    INSERT OR REPLACE INTO gpkg_geometry_columns 
+    (table_name, column_name, geometry_type_name, srs_id, z, m)
+    VALUES (?, 'geom', ?, ?, ?, 0)
+  """, [layer, geom_type, srs_id, z_flag])
+  
+  # Insert data row by row without prepared statements
+  if isnothing(table)
+    # Handle case with no attribute columns - just geometry
+    for (i, geom) in enumerate(geoms)
+      gpb_blob = create_gpb(geom, srs_id)
+      insert_sql = "INSERT INTO \"$layer\" (geom) VALUES (?)"
+      DBInterface.execute(db, insert_sql, [gpb_blob])
+    end
+  else
+    # Handle case with attribute columns
+    for (i, row) in enumerate(Tables.rows(table))
+      geom = geoms[i]
+      gpb_blob = create_gpb(geom, srs_id)
+      
+      # Build insert SQL dynamically
+      col_names = ["geom"; [string(col) for col in cols]]
+      col_placeholders = repeat(["?"], length(col_names))
+      insert_sql = "INSERT INTO \"$layer\" ($(join(col_names, ", "))) VALUES ($(join(col_placeholders, ", ")))"
+      
+      values = Any[gpb_blob]
+      for col in cols
+        val = getproperty(row, col)
+        push!(values, val)
+      end
+      
+      DBInterface.execute(db, insert_sql, values)
+    end
+  end
+  close(db)
+  
+  return fname
+end
\ No newline at end of file
diff --git a/src/precompile.jl b/src/precompile.jl
index 07a00a2..30968a9 100644
--- a/src/precompile.jl
+++ b/src/precompile.jl
@@ -12,9 +12,9 @@ using PrecompileTools
     load(joinpath(datadir, "points.parquet"))
     load(joinpath(datadir, "lines.parquet"))
     load(joinpath(datadir, "polygons.parquet"))
-    load(joinpath(datadir, "field.kml"))
     load(joinpath(datadir, "points.gpkg"))
     load(joinpath(datadir, "lines.gpkg"))
     load(joinpath(datadir, "polygons.gpkg"))
+    load(joinpath(datadir, "field.kml"))
   end
 end