Skip to content

Export API

Export modules — GeoTIFF, CSV, GeoJSON, QGIS project generation.

export_cog(data, output_path, **kwargs)

Export as Cloud-Optimized GeoTIFF (COG).

COGs enable efficient web streaming — stakeholders can explore visualizations without GIS software.

Parameters:

Name Type Description Default
data ndarray

2D or 3D array

 required
output_path str | Path

Output .tif path

 required
**kwargs Any

Passed to export_geotiff

 {}

Returns:

Type Description
str

Path to the written COG
Source code in dig/export/geotiff.py 
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def export_cog(
    data: np.ndarray,
    output_path: str | Path,
    **kwargs: Any,
) -> str:
    """Export as Cloud-Optimized GeoTIFF (COG).

    COGs enable efficient web streaming — stakeholders can explore
    visualizations without GIS software.

    Args:
        data: 2D or 3D array
        output_path: Output .tif path
        **kwargs: Passed to export_geotiff

    Returns:
        Path to the written COG
    """
    path = export_geotiff(data, output_path, **kwargs)

    # Re-open and convert to COG
    with rasterio.open(path) as src:
        profile = src.profile.copy()
        profile["driver"] = "COG"
        profile["compress"] = "DEFLATE"

        cog_path = Path(output_path).with_suffix(".cog.tif")
        with rasterio.open(cog_path, "w", **profile) as dst:
            for i in range(1, src.count + 1):
                dst.write(src.read(i), i)

    return str(cog_path)

export_csv(data, output_path, x_coords=None, y_coords=None, depth_labels=None)

Export processed data as CSV.

Parameters:

Name Type Description Default
data ndarray

2D (rows, cols) or 3D (bands, rows, cols) array

 required
output_path str | Path

Output .csv path

 required
x_coords ndarray | None

X coordinates for each column

 None
y_coords ndarray | None

Y coordinates for each row

 None
depth_labels list[str] | None

Labels for each band/depth

 None

Returns:

Type Description
str

Path to the written CSV
Source code in dig/export/csv_export.py 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
def export_csv(
    data: np.ndarray,
    output_path: str | Path,
    x_coords: np.ndarray | None = None,
    y_coords: np.ndarray | None = None,
    depth_labels: list[str] | None = None,
) -> str:
    """Export processed data as CSV.

    Args:
        data: 2D (rows, cols) or 3D (bands, rows, cols) array
        output_path: Output .csv path
        x_coords: X coordinates for each column
        y_coords: Y coordinates for each row
        depth_labels: Labels for each band/depth

    Returns:
        Path to the written CSV
    """
    output_path = Path(output_path)

    if data.ndim == 3:
        # Multi-band: one column per band
        n_bands, n_rows, n_cols = data.shape
        rows_list = []
        for i in range(n_rows):
            for j in range(n_cols):
                row = {
                    "row": i,
                    "col": j,
                }
                if x_coords is not None:
                    row["x"] = x_coords[j] if j < len(x_coords) else j
                if y_coords is not None:
                    row["y"] = y_coords[i] if i < len(y_coords) else i
                for b in range(n_bands):
                    label = (
                        depth_labels[b] if depth_labels and b < len(depth_labels) else f"band_{b}"
                    )
                    row[label] = data[b, i, j]
                rows_list.append(row)
        df = pd.DataFrame(rows_list)
    else:
        # 2D: simple grid
        df = pd.DataFrame(data)

    df.to_csv(output_path, index=False)
    return str(output_path)

export_geojson(features, output_path, crs_epsg=4326)

Export anomaly detections or survey boundaries as GeoJSON.

Parameters:

Name Type Description Default
features list[dict]

List of GeoJSON Feature objects

 required
output_path str | Path

Output .geojson path

 required
crs_epsg int

EPSG code

 4326

Returns:

Type Description
str

Path to the written GeoJSON
Source code in dig/export/geojson.py 
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def export_geojson(
    features: list[dict],
    output_path: str | Path,
    crs_epsg: int = 4326,
) -> str:
    """Export anomaly detections or survey boundaries as GeoJSON.

    Args:
        features: List of GeoJSON Feature objects
        output_path: Output .geojson path
        crs_epsg: EPSG code

    Returns:
        Path to the written GeoJSON
    """
    output_path = Path(output_path)

    geojson = {
        "type": "FeatureCollection",
        "features": features,
        "crs": {
            "type": "name",
            "properties": {"name": f"EPSG:{crs_epsg}"},
        },
    }

    with open(output_path, "w") as f:
        json.dump(geojson, f, indent=2)

    return str(output_path)

export_geotiff(data, output_path, crs_epsg=4326, origin_easting=0.0, origin_northing=0.0, pixel_size=1.0, rotation_deg=0.0, band_descriptions=None, nodata=None)

Export processed geophysical data as a georeferenced GeoTIFF.

Supports multi-band export (each band = one depth slice).

Parameters:

Name Type Description Default
data ndarray

2D (rows, cols) or 3D (bands, rows, cols) array

 required
output_path str | Path

Output .tif path

 required
crs_epsg int

EPSG code for the CRS

 4326
origin_easting float

Easting of grid origin

 0.0
origin_northing float

Northing of grid origin

 0.0
pixel_size float

Pixel resolution in CRS units

 1.0
rotation_deg float

Clockwise rotation from north

 0.0
band_descriptions Optional[list[str]]

Per-band description strings

 None
nodata Optional[float]

NoData value (auto-detected if None)

 None

Returns:

Type Description
str

Path to the written GeoTIFF
Source code in dig/export/geotiff.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
def export_geotiff(
    data: np.ndarray,
    output_path: str | Path,
    crs_epsg: int = 4326,
    origin_easting: float = 0.0,
    origin_northing: float = 0.0,
    pixel_size: float = 1.0,
    rotation_deg: float = 0.0,
    band_descriptions: Optional[list[str]] = None,
    nodata: Optional[float] = None,
) -> str:
    """Export processed geophysical data as a georeferenced GeoTIFF.

    Supports multi-band export (each band = one depth slice).

    Args:
        data: 2D (rows, cols) or 3D (bands, rows, cols) array
        output_path: Output .tif path
        crs_epsg: EPSG code for the CRS
        origin_easting: Easting of grid origin
        origin_northing: Northing of grid origin
        pixel_size: Pixel resolution in CRS units
        rotation_deg: Clockwise rotation from north
        band_descriptions: Per-band description strings
        nodata: NoData value (auto-detected if None)

    Returns:
        Path to the written GeoTIFF
    """
    output_path = Path(output_path)

    # Handle 2D → 3D
    if data.ndim == 2:
        data = data[np.newaxis, :, :]

    n_bands, n_rows, n_cols = data.shape

    # Build affine transform
    transform = _build_affine(
        origin_easting,
        origin_northing,
        pixel_size,
        rotation_deg,
        n_cols,
        n_rows,
    )

    # CRS
    crs = CRS.from_epsg(crs_epsg)

    # Auto-detect nodata
    if nodata is None:
        if np.issubdtype(data.dtype, np.floating):
            nodata = np.nan
        else:
            nodata = -9999

    profile = {
        "driver": "GTiff",
        "height": n_rows,
        "width": n_cols,
        "count": n_bands,
        "dtype": data.dtype,
        "crs": crs,
        "transform": transform,
        "nodata": nodata,
        "compress": "lzw",
        "tiled": True,
        "blockxsize": 256,
        "blockysize": 256,
    }

    with rasterio.open(output_path, "w", **profile) as dst:
        for i in range(n_bands):
            band_data = data[i]
            if nodata is not None and np.isnan(nodata):
                band_data = np.where(np.isnan(band_data), nodata, band_data)
            dst.write(band_data, i + 1)

            # Band description
            if band_descriptions and i < len(band_descriptions):
                dst.set_band_description(i + 1, band_descriptions[i])

    return str(output_path)

generate_qgz(raster_paths, output_path, project_title='DIG Geophysical Survey', crs_epsg=4326)

Generate a QGIS project file (.qgs XML).

Configures a complete QGIS 3.x project with layer symbology, global extent, and proper project tree structuring.

Parameters:

Name Type Description Default
raster_paths list[str]

Paths to GeoTIFF files to include

 required
output_path str | Path

Output .qgs path

 required
project_title str

Project title

 'DIG Geophysical Survey'
crs_epsg int

EPSG code for the project CRS

 4326

Returns:

Type Description
str

Path to the written .qgs file
Source code in dig/export/qgis.py 
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
def generate_qgz(
    raster_paths: list[str],
    output_path: str | Path,
    project_title: str = "DIG Geophysical Survey",
    crs_epsg: int = 4326,
) -> str:
    """Generate a QGIS project file (.qgs XML).

    Configures a complete QGIS 3.x project with layer symbology, global extent,
    and proper project tree structuring.

    Args:
        raster_paths: Paths to GeoTIFF files to include
        output_path: Output .qgs path
        project_title: Project title
        crs_epsg: EPSG code for the project CRS

    Returns:
        Path to the written .qgs file
    """
    output_path = Path(output_path)

    xmin, ymin, xmax, ymax = float("inf"), float("inf"), float("-inf"), float("-inf")

    projectlayers_xml = ""
    layertree_xml = ""

    for i, rpath_str in enumerate(raster_paths):
        rpath = Path(rpath_str)
        layer_id = f"dig_layer_{i}"
        layer_name = rpath.stem

        vmin, vmax = 0.0, 255.0

        try:
            with rasterio.open(rpath) as src:
                bounds = src.bounds
                xmin = min(xmin, bounds.left)
                ymin = min(ymin, bounds.bottom)
                xmax = max(xmax, bounds.right)
                ymax = max(ymax, bounds.top)

                band = src.read(1)
                nodata = src.nodata

                if nodata is not None:
                    if np.issubdtype(band.dtype, np.floating) and np.isnan(nodata):
                        valid = band[~np.isnan(band)]
                    else:
                        valid = band[band != nodata]
                else:
                    if np.issubdtype(band.dtype, np.floating):
                        valid = band[~np.isnan(band)]
                    else:
                        valid = band

                if valid.size > 0:
                    vmin = float(np.percentile(valid, 2))
                    vmax = float(np.percentile(valid, 98))
        except Exception:
            pass

        projectlayers_xml += f"""
    <maplayer type="raster" hasScaleBasedVisibilityFlag="0" maxScale="0" minScale="1e+08" autoRefreshTime="0" autoRefreshEnabled="0" refreshOnNotifyEnabled="0" refreshOnNotifyMessage="">
      <id>{layer_id}</id>
      <datasource>{rpath.resolve()}</datasource>
      <layername>{layer_name}</layername>
      <srs>
        <spatialrefsys>
          <authid>EPSG:{crs_epsg}</authid>
        </spatialrefsys>
      </srs>
      <provider>gdal</provider>
      <pipe>
        <rasterrenderer type="singlebandpseudocolor" band="1" alphaBand="-1" classificationMin="{vmin}" classificationMax="{vmax}" opacity="1">
          <rasterTransparency/>
          <minMaxOrigin>
            <limits>MinMax</limits>
            <extent>WholeRaster</extent>
            <statAccuracy>Estimated</statAccuracy>
            <cumulativeCutLower>0.02</cumulativeCutLower>
            <cumulativeCutUpper>0.98</cumulativeCutUpper>
            <stdDevFactor>2</stdDevFactor>
          </minMaxOrigin>
          <rastershader>
            <colorrampshader colorRampType="INTERPOLATED" clip="0" classificationMode="1" maximumValue="{vmax}" minimumValue="{vmin}">
              <colorramp type="gradient" name="[source]">
                <prop k="color1" v="43,131,186,255"/>
                <prop k="color2" v="215,25,28,255"/>
                <prop k="discrete" v="0"/>
                <prop k="rampType" v="gradient"/>
              </colorramp>
              <item alpha="255" value="{vmin}" label="{vmin:.1f}" color="#2b83ba"/>
              <item alpha="255" value="{(vmin + vmax) / 2}" label="{(vmin + vmax) / 2:.1f}" color="#ffffbf"/>
              <item alpha="255" value="{vmax}" label="{vmax:.1f}" color="#d7191c"/>
            </colorrampshader>
          </rastershader>
        </rasterrenderer>
      </pipe>
    </maplayer>"""

        layertree_xml += f"""
    <layer-tree-layer name="{layer_name}" providerKey="gdal" source="{rpath.resolve()}" checked="Qt::Checked" id="{layer_id}" expanded="1">
      <customproperties/>
    </layer-tree-layer>"""

    if math.isinf(xmin):
        xmin, ymin, xmax, ymax = 0, 0, 1000, 1000

    qgs_content = f"""<!DOCTYPE qgis PUBLIC 'http://mrcc.com/qgis.dtd' 'SYSTEM'>
<qgis version="3.22.0" projectname="{project_title}">
  <title>{project_title}</title>
  <projectlayers>{projectlayers_xml}
  </projectlayers>
  <layer-tree-group>
    <customproperties/>{layertree_xml}
  </layer-tree-group>
  <mapcanvas>
    <units>meters</units>
    <extent>
      <xmin>{xmin}</xmin>
      <ymin>{ymin}</ymin>
      <xmax>{xmax}</xmax>
      <ymax>{ymax}</ymax>
    </extent>
    <projections>
      <crs>
        <spatialrefsys>
          <authid>EPSG:{crs_epsg}</authid>
        </spatialrefsys>
      </crs>
    </projections>
  </mapcanvas>
</qgis>"""

    output_path.write_text(qgs_content)
    return str(output_path)