Skip to content

Raster to DGGS

Raster to DGGS conversion functions.

This submodule provides functions to convert raster data to various discrete global grid systems (DGGS).

raster2olc_cli()

Command line interface for raster to OLC conversion

Source code in vgrid/conversion/raster2dggs/raster2olc.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def raster2olc_cli():
    """Command line interface for raster to OLC conversion"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to OLC/ Google Plus Code DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")

    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        choices=olc_res,
        default=None,
        help="OLC resolution",
    )

    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )

    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format

    if not os.path.exists(raster):
        print(f"Error: The file {raster} does not exist.")
        return

    result = raster2olc(raster, resolution, output_format)
    if output_format in STRUCTURED_FORMATS:
        print(result)

raster2qtm_cli()

Command line interface for raster2qtm conversion

Source code in vgrid/conversion/raster2dggs/raster2qtm.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def raster2qtm_cli():
    """Command line interface for raster2qtm conversion"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to QTM DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")
    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        default=None,
        help=f"QTM resolution [{min_res}..{max_res}]",
    )
    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )

    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format

    if not os.path.exists(raster):
        print(f"Error: The file {raster} does not exist.")
        return

    result = raster2qtm(raster, resolution, output_format)
    if output_format in STRUCTURED_FORMATS:
        print(result)

raster2rhealpix_cli()

Command line interface for raster2rhealpix

Source code in vgrid/conversion/raster2dggs/raster2rhealpix.py 
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
def raster2rhealpix_cli():
    """Command line interface for raster2rhealpix"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to rHEALPix DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")
    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        default=None,
        help=f"rHEALPix resolution [{min_res}..{max_res}]. Required when topology=False, auto-calculated when topology=True",
    )
    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )
    parser.add_argument(
        "-fix",
        "--fix_antimeridian",
        type=str,
        choices=[
            "shift",
            "shift_balanced",
            "shift_west",
            "shift_east",
            "split",
            "none",
        ],
        default=None,
        help="Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none",
    )
    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format
    fix_antimeridian = args.fix_antimeridian
    if not os.path.exists(raster):
        raise FileNotFoundError(f"The file {raster} does not exist.")

    result = raster2rhealpix(
        raster, resolution, output_format, fix_antimeridian=fix_antimeridian
    )
    if output_format in STRUCTURED_FORMATS:
        print(result)

Raster to H3 Module

This module provides functionality to convert raster data to H3 (Hierarchical Hexagonal Grid) DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_h3_resolution(raster_path)

Automatically determine the optimal H3 resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate H3 resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal H3 resolution level

Examples

cell_size, resolution = get_nearest_h3_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2h3.py 
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
def get_nearest_h3_resolution(raster_path):
    """
    Automatically determine the optimal H3 resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate H3 resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal H3 resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_h3_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    # Check resolutions from 0 to 15
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        avg_area = h3.average_hexagon_area(res, unit="m^2")
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2h3(raster_path, resolution=None, output_format='gpd', fix_antimeridian=None)

Convert raster data to H3 DGGS format.

Parameters:

Name Type Description Default
raster_path str

Path to the raster file

 required
resolution int

H3 resolution [0..15]. If None, automatically determined

 None
output_format str

Output format. Options: - None: Returns GeoPandas GeoDataFrame (default) - "gpd": Returns GeoPandas GeoDataFrame - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

 'gpd'
Source code in vgrid/conversion/raster2dggs/raster2h3.py 
 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
155
156
157
158
159
160
161
162
163
164
165
def raster2h3(raster_path, resolution=None, output_format="gpd", fix_antimeridian=None):
    """
    Convert raster data to H3 DGGS format.

    Args:
        raster_path (str): Path to the raster file
        resolution (int, optional): H3 resolution [0..15]. If None, automatically determined
        output_format (str, optional): Output format. Options:
            - None: Returns GeoPandas GeoDataFrame (default)
            - "gpd": Returns GeoPandas GeoDataFrame
            - "csv": Returns CSV file path
            - "geojson": Returns GeoJSON file path
            - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
            - "parquet": Returns Parquet file path
            - "shapefile"/"shp": Returns Shapefile file path
            - "gpkg"/"geopackage": Returns GeoPackage file path
    Returns:
        Various formats based on output_format parameter
    Raises:
        ValueError: If resolution is not in valid range [0..15]
        ImportError: If required dependencies are not available for specific formats
    """
    # Step 1: Determine the nearest H3 resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_h3_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest H3 resolution determined: {resolution}")
    else:
        resolution = validate_h3_resolution(resolution)

    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per H3 cell
    h3_ids_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            h3_id = h3.latlng_to_cell(lat, lon, resolution)
            if h3_id not in h3_ids_band_values:
                vals = raster_data[:, int(row), int(col)]
                # Convert NumPy scalars to native Python types
                h3_ids_band_values[h3_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for h3_id, band_values in tqdm(
        h3_ids_band_values.items(), desc="Converting raster to H3", unit=" cells"
    ):
        cell_polygon = h32geo(h3_id, fix_antimeridian=fix_antimeridian)
        num_edges = 6
        if h3.is_pentagon(h3_id):
            num_edges = 5
        base_props = geodesic_dggs_to_geoseries(
            "h3", h3_id, resolution, cell_polygon, num_edges
        )
        band_properties = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_properties)
        properties.append(base_props)

    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2h3" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to S2 Module

This module provides functionality to convert raster data to S2 DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_s2_resolution(raster_path)

Automatically determine the optimal S2 resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate S2 resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal S2 resolution level

Examples

cell_size, resolution = get_nearest_s2_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2s2.py 
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
def get_nearest_s2_resolution(raster_path):
    """
    Automatically determine the optimal S2 resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate S2 resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal S2 resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_s2_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    # Check resolutions from 0 to 15
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = s2_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2s2(raster_path, resolution=None, output_format='gpd', fix_antimeridian=None)

Convert raster data to S2 DGGS format.

Converts raster data to S2 DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to an S2 cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional S2 resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-30. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path fix_antimeridian : str, optional Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none Defaults to None when omitted. Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents an S2 cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2s2("data.tif") print(f"Converted {len(result)} S2 cells")

Convert with specific resolution

result = raster2s2("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2s2("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2s2.py 
 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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
def raster2s2(raster_path, resolution=None, output_format="gpd", fix_antimeridian=None):
    """
    Convert raster data to S2 DGGS format.

    Converts raster data to S2 DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to an S2 cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        S2 resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-30.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    fix_antimeridian : str, optional
        Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none
        Defaults to None when omitted.
    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents an S2 cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2s2("data.tif")
    >>> print(f"Converted {len(result)} S2 cells")

    >>> # Convert with specific resolution
    >>> result = raster2s2("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2s2("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest s2 resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_s2_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest S2 resolution determined: {resolution}")
    else:
        resolution = validate_s2_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per S2 cell
    s2_tokens_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            lat_lng = s2.LatLng.from_degrees(lat, lon)
            s2_id = s2.CellId.from_lat_lng(lat_lng)
            s2_id = s2_id.parent(resolution)
            s2_token = s2.CellId.to_token(s2_id)
            if s2_token not in s2_tokens_band_values:
                vals = raster_data[:, int(row), int(col)]
                s2_tokens_band_values[s2_token] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    # Build GeoDataFrame as the base
    properties = []
    for s2_token, band_values in tqdm(
        s2_tokens_band_values.items(), desc="Converting raster to S2", unit=" cells"
    ):
        cell_polygon = s22geo(s2_token, fix_antimeridian=fix_antimeridian)
        num_edges = 4
        centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
            geodesic_dggs_metrics(cell_polygon, num_edges)
        )
        base_props = {
            "s2": s2_token,
            "resolution": resolution,
            "center_lat": centroid_lat,
            "center_lon": centroid_lon,
            "avg_edge_len": avg_edge_len,
            "cell_area": cell_area,
            "cell_perimeter": cell_perimeter,
            "geometry": cell_polygon,
        }
        band_properties = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_properties)
        properties.append(base_props)
    # Build GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2s2" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to A5 Module

This module provides functionality to convert raster data to A5 DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_a5_resolution(raster_path)

Automatically determine the optimal A5 resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate A5 resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal A5 resolution level

Examples

cell_size, resolution = get_nearest_a5_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2a5.py 
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
def get_nearest_a5_resolution(raster_path):
    """
    Automatically determine the optimal A5 resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate A5 resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal A5 resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_a5_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res
    for res in range(min_res, max_res + 1):
        avg_area = cell_area(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2a5(raster_path, resolution=None, output_format='gpd', options=None, split_antimeridian=False)

Convert raster data to A5 DGGS format.

Converts raster data to A5 (Adaptive 5) DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to an A5 cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional A5 resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-15. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path options : dict, optional Options for a52geo. split_antimeridian : bool, optional When True, apply antimeridian fixing to the resulting polygons. Defaults to False when None or omitted. Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents an A5 cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2a5("data.tif") print(f"Converted {len(result)} A5 cells")

Convert with specific resolution

result = raster2a5("data.tif", resolution=5)

Convert to GeoJSON file

result = raster2a5("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2a5.py 
 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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
def raster2a5(
    raster_path,
    resolution=None,
    output_format="gpd",
    options=None,
    split_antimeridian=False,
):
    """
    Convert raster data to A5 DGGS format.

    Converts raster data to A5 (Adaptive 5) DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to an A5 cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        A5 resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-15.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    options : dict, optional
        Options for a52geo.
    split_antimeridian : bool, optional
        When True, apply antimeridian fixing to the resulting polygons.
        Defaults to False when None or omitted.
    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents an A5 cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2a5("data.tif")
    >>> print(f"Converted {len(result)} A5 cells")

    >>> # Convert with specific resolution
    >>> result = raster2a5("data.tif", resolution=5)

    >>> # Convert to GeoJSON file
    >>> result = raster2a5("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest a5 resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_a5_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest A5 resolution determined: {resolution}")
    else:
        resolution = validate_a5_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per A5 cell
    a5_hexs_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            a5_hex = latlon2a5(lat, lon, resolution)
            if a5_hex not in a5_hexs_band_values:
                vals = raster_data[:, int(row), int(col)]
                a5_hexs_band_values[a5_hex] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    # Build GeoDataFrame as the base
    properties = []
    for a5_hex, band_values in tqdm(
        a5_hexs_band_values.items(), desc="Converting raster to A5", unit=" cells"
    ):
        try:
            cell_polygon = a52geo(
                a5_hex, options, split_antimeridian=split_antimeridian
            )
            num_edges = 5
            centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
                geodesic_dggs_metrics(cell_polygon, num_edges)
            )
            base_props = {
                "a5": a5_hex,
                "resolution": resolution,
                "center_lat": centroid_lat,
                "center_lon": centroid_lon,
                "avg_edge_len": avg_edge_len,
                "cell_area": cell_area,
                "cell_perimeter": cell_perimeter,
                "geometry": cell_polygon,
            }
            band_properties = {
                f"band_{i + 1}": band_values[i] for i in range(band_count)
            }
            base_props.update(band_properties)
            properties.append(base_props)
        except Exception:
            continue

    # Build GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2a5" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to RHEALPix Module

This module provides functionality to convert raster data to RHEALPix (Rectified HEALPix) DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_rhealpix_resolution(raster_path)

Automatically determine the optimal RHEALPix resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate RHEALPix resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal RHEALPix resolution level

Examples

cell_size, resolution = get_nearest_rhealpix_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2rhealpix.py 
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
def get_nearest_rhealpix_resolution(raster_path):
    """
    Automatically determine the optimal RHEALPix resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate RHEALPix resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal RHEALPix resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_rhealpix_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = rhealpix_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2rhealpix(raster_path, resolution=None, output_format='gpd', fix_antimeridian=None)

Convert raster data to RHEALPix DGGS format.

Converts raster data to RHEALPix (Rectified HEALPix) DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a RHEALPix cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional RHEALPix resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-30. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path fix_antimeridian : Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none When True, apply antimeridian fixing to the resulting polygons. Defaults to False when None or omitted.

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a RHEALPix cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2rhealpix("data.tif") print(f"Converted {len(result)} RHEALPix cells")

Convert with specific resolution

result = raster2rhealpix("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2rhealpix("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2rhealpix.py 
 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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
def raster2rhealpix(
    raster_path, resolution=None, output_format="gpd", fix_antimeridian=None
):
    """
    Convert raster data to RHEALPix DGGS format.

    Converts raster data to RHEALPix (Rectified HEALPix) DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a RHEALPix cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        RHEALPix resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-30.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    fix_antimeridian : Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none
        When True, apply antimeridian fixing to the resulting polygons.
        Defaults to False when None or omitted.

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a RHEALPix cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2rhealpix("data.tif")
    >>> print(f"Converted {len(result)} RHEALPix cells")

    >>> # Convert with specific resolution
    >>> result = raster2rhealpix("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2rhealpix("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest rhealpix resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_rhealpix_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest rHEALPix resolution determined: {resolution}")
    else:
        resolution = validate_rhealpix_resolution(resolution)

    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per rHEALPix cell
    rhealpix_ids_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            point = (lon, lat)
            rhealpix_cell = rhealpix_dggs.cell_from_point(
                resolution, point, plane=False
            )
            rhealpix_id = str(rhealpix_cell)
            if rhealpix_id not in rhealpix_ids_band_values:
                vals = raster_data[:, int(row), int(col)]
                rhealpix_ids_band_values[rhealpix_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    # Build GeoDataFrame as the base
    properties = []
    for rhealpix_id, band_values in tqdm(
        rhealpix_ids_band_values.items(),
        desc="Converting raster to rHEALPix",
        unit=" cells",
    ):
        cell_polygon = rhealpix2geo(rhealpix_id, fix_antimeridian=fix_antimeridian)
        rhealpix_uids = (rhealpix_id[0],) + tuple(map(int, rhealpix_id[1:]))
        rhealpix_cell = rhealpix_dggs.cell(rhealpix_uids)
        num_edges = 4
        if rhealpix_cell.ellipsoidal_shape() == "dart":
            num_edges = 3
        centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
            geodesic_dggs_metrics(cell_polygon, num_edges)
        )
        base_props = {
            "rhealpix": rhealpix_id,
            "resolution": resolution,
            "center_lat": centroid_lat,
            "center_lon": centroid_lon,
            "avg_edge_len": avg_edge_len,
            "cell_area": cell_area,
            "cell_perimeter": cell_perimeter,
            "geometry": cell_polygon,
        }
        band_properties = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_properties)
        properties.append(base_props)
    # Build GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2rhealpix" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

raster2rhealpix_cli()

Command line interface for raster2rhealpix

Source code in vgrid/conversion/raster2dggs/raster2rhealpix.py 
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
def raster2rhealpix_cli():
    """Command line interface for raster2rhealpix"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to rHEALPix DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")
    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        default=None,
        help=f"rHEALPix resolution [{min_res}..{max_res}]. Required when topology=False, auto-calculated when topology=True",
    )
    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )
    parser.add_argument(
        "-fix",
        "--fix_antimeridian",
        type=str,
        choices=[
            "shift",
            "shift_balanced",
            "shift_west",
            "shift_east",
            "split",
            "none",
        ],
        default=None,
        help="Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none",
    )
    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format
    fix_antimeridian = args.fix_antimeridian
    if not os.path.exists(raster):
        raise FileNotFoundError(f"The file {raster} does not exist.")

    result = raster2rhealpix(
        raster, resolution, output_format, fix_antimeridian=fix_antimeridian
    )
    if output_format in STRUCTURED_FORMATS:
        print(result)

Raster to ISEA4T Module

This module provides functionality to convert raster data to ISEA4T DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_isea4t_resolution(raster_path)

Automatically determine the optimal ISEA4T resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate ISEA4T resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal ISEA4T resolution level

Examples

cell_size, resolution = get_nearest_isea4t_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2isea4t.py 
 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
def get_nearest_isea4t_resolution(raster_path):
    """
    Automatically determine the optimal ISEA4T resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate ISEA4T resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal ISEA4T resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_isea4t_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = isea4t_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2isea4t(raster_path, resolution=None, output_format='gpd', fix_antimeridian=None)

Convert raster data to ISEA4T DGGS format.

Converts raster data to ISEA4T DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to an ISEA4T cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional ISEA4T resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-39. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path fix_antimeridian : str, optional Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents an ISEA4T cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2isea4t("data.tif") print(f"Converted {len(result)} ISEA4T cells")

Convert with specific resolution

result = raster2isea4t("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2isea4t("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2isea4t.py 
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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
def raster2isea4t(
    raster_path, resolution=None, output_format="gpd", fix_antimeridian=None
):
    """
    Convert raster data to ISEA4T DGGS format.

    Converts raster data to ISEA4T DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to an ISEA4T cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        ISEA4T resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-39.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    fix_antimeridian : str, optional
        Antimeridian fixing method: shift, shift_balanced, shift_west, shift_east, split, none

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents an ISEA4T cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2isea4t("data.tif")
    >>> print(f"Converted {len(result)} ISEA4T cells")

    >>> # Convert with specific resolution
    >>> result = raster2isea4t("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2isea4t("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest isea4t resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_isea4t_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest ISEA4T resolution determined: {resolution}")
    else:
        resolution = validate_isea4t_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per ISEA4T cell
    isea4t_ids_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            max_accuracy = ISEA4T_RES_ACCURACY_DICT[39]
            lat_long_point = LatLongPoint(lat, lon, max_accuracy)
            isea4t_cell_max_accuracy = isea4t_dggs.convert_point_to_dggs_cell(
                lat_long_point
            )
            cell_id_len = resolution + 2
            isea4t_cell = DggsCell(isea4t_cell_max_accuracy._cell_id[:cell_id_len])
            isea4t_id = isea4t_cell._cell_id
            if isea4t_id not in isea4t_ids_band_values:
                vals = raster_data[:, int(row), int(col)]
                isea4t_ids_band_values[isea4t_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for isea4t_id, band_values in tqdm(
        isea4t_ids_band_values.items(),
        desc="Converting raster to ISEA4T",
        unit=" cells",
    ):
        cell_polygon = isea4t2geo(isea4t_id, fix_antimeridian=fix_antimeridian)
        num_edges = 3
        centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
            geodesic_dggs_metrics(cell_polygon, num_edges)
        )
        base_props = {
            "isea4t": isea4t_id,
            "resolution": resolution,
            "center_lat": centroid_lat,
            "center_lon": centroid_lon,
            "avg_edge_len": avg_edge_len,
            "cell_area": cell_area,
            "cell_perimeter": cell_perimeter,
            "geometry": cell_polygon,
        }
        band_properties = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_properties)
        properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2isea4t" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to DGGAL Module

This module provides functionality to convert raster data to DGGAL (Discrete Global Grids with Adaptive Localization) DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_dggal_resolution(dggs_type, raster_path)

Automatically determine the optimal DGGAL resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate DGGAL resolution that best matches the raster's spatial resolution for the specified DGGS type.

Parameters

dggs_type : str DGGAL DGGS type (e.g., "isea3h", "isea4t", "rhealpix"). raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal DGGAL resolution level

Examples

cell_size, resolution = get_nearest_dggal_resolution("isea3h", "data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2dggal.py 
 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
def get_nearest_dggal_resolution(dggs_type, raster_path):
    """
    Automatically determine the optimal DGGAL resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate DGGAL resolution
    that best matches the raster's spatial resolution for the specified DGGS type.

    Parameters
    ----------
    dggs_type : str
        DGGAL DGGS type (e.g., "isea3h", "isea4t", "rhealpix").
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal DGGAL resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_dggal_resolution("isea3h", "data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * math.cos(
                math.radians(center_latitude)
            )

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    min_res = int(DGGAL_TYPES[dggs_type]["min_res"])
    max_res = int(DGGAL_TYPES[dggs_type]["max_res"])
    nearest_resolution = min_res
    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = dggal_metrics(dggs_type, res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = math.fabs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2dggal(dggs_type, raster_path, resolution=None, output_format='gpd', split_antimeridian=False)

Convert raster data to DGGAL DGGS format.

Converts raster data to DGGAL DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a DGGAL cell and the first sample value per cell is preserved.

Parameters

dggs_type : str DGGAL DGGS type (e.g., "isea3h", "isea4t", "rhealpix"). raster_path : str Path to the raster file to convert. resolution : int, optional DGGAL resolution level. If None, automatically determined based on raster pixel size. Valid range depends on the DGGS type. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path split_antimeridian : bool, optional When True, apply antimeridian fixing to the resulting polygons. Defaults to False when None or omitted.

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a DGGAL cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2dggal("isea3h", "data.tif") print(f"Converted {len(result)} DGGAL cells")

Convert with specific resolution

result = raster2dggal("isea3h", "data.tif", resolution=5)

Convert to GeoJSON file

result = raster2dggal("isea3h", "data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2dggal.py 
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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def raster2dggal(
    dggs_type: str,
    raster_path,
    resolution: int | None = None,
    output_format: str = "gpd",
    split_antimeridian: bool = False,
):
    """
    Convert raster data to DGGAL DGGS format.

    Converts raster data to DGGAL DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a DGGAL cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    dggs_type : str
        DGGAL DGGS type (e.g., "isea3h", "isea4t", "rhealpix").
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        DGGAL resolution level. If None, automatically determined based on raster pixel size.
        Valid range depends on the DGGS type.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    split_antimeridian : bool, optional
        When True, apply antimeridian fixing to the resulting polygons.
        Defaults to False when None or omitted.

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a DGGAL cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2dggal("isea3h", "data.tif")
    >>> print(f"Converted {len(result)} DGGAL cells")

    >>> # Convert with specific resolution
    >>> result = raster2dggal("isea3h", "data.tif", resolution=5)

    >>> # Convert to GeoJSON file
    >>> result = raster2dggal("isea3h", "data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    dggs_type = validate_dggal_type(dggs_type)
    # Auto-select resolution if not provided
    if resolution is None:
        cell_size, resolution = get_nearest_dggal_resolution(dggs_type, raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest {dggs_type.upper()} resolution determined: {resolution}")
    else:
        resolution = validate_dggal_resolution(dggs_type, resolution)

    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per DGGAL cell
    zone_ids_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            try:
                zone_id = latlon2dggal(dggs_type, lat, lon, resolution)
                if zone_id not in zone_ids_band_values:
                    vals = raster_data[:, int(row), int(col)]
                    zone_ids_band_values[zone_id] = [
                        (v.item() if hasattr(v, "item") else v) for v in vals
                    ]
            except Exception:
                continue
    # Build GeoDataFrame as the base
    properties = []
    for zone_id, band_values in tqdm(
        zone_ids_band_values.items(), desc="Converting raster to DGGAL", unit=" cells"
    ):
        try:
            # Get zone object to get resolution and edge count
            dggs_class_name = DGGAL_TYPES[dggs_type]["class_name"]
            dggrs = globals()[dggs_class_name]()
            zone = dggrs.getZoneFromTextID(zone_id)
            num_edges = dggrs.countZoneEdges(zone)

            # Convert zone to geometry using dggal2geo
            cell_polygon = dggal2geo(
                dggs_type, zone_id, split_antimeridian=split_antimeridian
            )

            centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
                geodesic_dggs_metrics(cell_polygon, num_edges)
            )
            base_props = {
                f"dggal_{dggs_type}": zone_id,
                "resolution": resolution,
                "center_lat": centroid_lat,
                "center_lon": centroid_lon,
                "avg_edge_len": avg_edge_len,
                "cell_area": cell_area,
                "cell_perimeter": cell_perimeter,
                "geometry": cell_polygon,
            }
            band_properties = {
                f"band_{i + 1}": band_values[i] for i in range(band_count)
            }
            base_props.update(band_properties)
            properties.append(base_props)
        except Exception:
            continue

    # Build GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2dggal" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to DGGRID Module

This module provides functionality to convert raster data to DGGRID DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_dggrid_resolution(dggrid_instance, dggs_type, raster_path)

Automatically determine the optimal DGGRID resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate DGGRID resolution that best matches the raster's spatial resolution for the specified DGGS type.

Parameters

dggrid_instance : object DGGRID instance for processing. dggs_type : str DGGRID DGGS type (e.g., "ISEA7H", "ISEA4T"). raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal DGGRID resolution level

Examples

cell_size, resolution = get_nearest_dggrid_resolution(instance, "ISEA7H", "data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2dggrid.py 
 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
def get_nearest_dggrid_resolution(dggrid_instance, dggs_type, raster_path):
    """
    Automatically determine the optimal DGGRID resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate DGGRID resolution
    that best matches the raster's spatial resolution for the specified DGGS type.

    Parameters
    ----------
    dggrid_instance : object
        DGGRID instance for processing.
    dggs_type : str
        DGGRID DGGS type (e.g., "ISEA7H", "ISEA4T").
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal DGGRID resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_dggrid_resolution(instance, "ISEA7H", "data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * math.cos(
                math.radians(center_latitude)
            )

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    min_res = int(DGGRID_TYPES[dggs_type]["min_res"])
    max_res = int(DGGRID_TYPES[dggs_type]["max_res"])
    nearest_resolution = min_res

    try:
        # Get stats with area in m^2 to compare directly with raster cell size
        grid_stats = dggridstats(dggrid_instance, dggs_type, unit="m")
        for res in range(min_res, max_res + 1):
            res_stats = grid_stats[grid_stats["resolution"] == res]
            if res_stats.empty:
                continue
            avg_area_m2 = res_stats["area_m2"].iloc[0]
            if avg_area_m2 < MIN_CELL_AREA:
                break
            diff = math.fabs(avg_area_m2 - cell_size)
            if diff < min_diff:
                min_diff = diff
                nearest_resolution = res
    except Exception:
        # Fallback to using min_res if grid stats fail
        nearest_resolution = min_res

    return cell_size, nearest_resolution

raster2dggrid(dggrid_instance, dggs_type, raster_path, resolution=None, output_format='gpd', split_antimeridian=False, aggregate=False)

Convert raster data to DGGRID DGGS format.

Converts raster data to DGGRID DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a DGGRID cell and the first sample value per cell is preserved.

Parameters

dggrid_instance : object DGGRID instance for processing. dggs_type : str DGGRID DGGS type (e.g., "ISEA7H", "ISEA4T"). raster_path : str Path to the raster file to convert. resolution : int, optional DGGRID resolution level. If None, automatically determined based on raster pixel size. Valid range depends on the DGGS type. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path split_antimeridian : bool, optional When True, apply antimeridian fixing to the resulting polygons. Defaults to False when None or omitted. aggregate : bool, optional When True, aggregate the resulting polygons. Defaults to False when None or omitted. Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a DGGRID cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2dggrid(instance, "ISEA7H", "data.tif") print(f"Converted {len(result)} DGGRID cells")

Convert with specific resolution

result = raster2dggrid(instance, "ISEA7H", "data.tif", resolution=5)

Convert to GeoJSON file

result = raster2dggrid(instance, "ISEA7H", "data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2dggrid.py 
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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
def raster2dggrid(
    dggrid_instance,
    dggs_type: str,
    raster_path,
    resolution: int | None = None,
    output_format: str = "gpd",
    split_antimeridian: bool = False,
    aggregate: bool = False,
):
    """
    Convert raster data to DGGRID DGGS format.

    Converts raster data to DGGRID DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a DGGRID cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    dggrid_instance : object
        DGGRID instance for processing.
    dggs_type : str
        DGGRID DGGS type (e.g., "ISEA7H", "ISEA4T").
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        DGGRID resolution level. If None, automatically determined based on raster pixel size.
        Valid range depends on the DGGS type.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path
    split_antimeridian : bool, optional
        When True, apply antimeridian fixing to the resulting polygons.
        Defaults to False when None or omitted.
    aggregate : bool, optional
        When True, aggregate the resulting polygons.
        Defaults to False when None or omitted.
    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a DGGRID cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2dggrid(instance, "ISEA7H", "data.tif")
    >>> print(f"Converted {len(result)} DGGRID cells")

    >>> # Convert with specific resolution
    >>> result = raster2dggrid(instance, "ISEA7H", "data.tif", resolution=5)

    >>> # Convert to GeoJSON file
    >>> result = raster2dggrid(instance, "ISEA7H", "data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    dggs_type = validate_dggrid_type(dggs_type)

    # Auto-select resolution if not provided
    if resolution is None:
        cell_size, resolution = get_nearest_dggrid_resolution(
            dggrid_instance, dggs_type, raster_path
        )
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest {dggs_type.upper()} resolution determined: {resolution}")
    else:
        resolution = validate_dggrid_resolution(dggs_type, resolution)

    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per DGGRID cell
    dggrid_ids_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            try:
                dggrid_id = latlon2dggrid(
                    dggrid_instance, dggs_type, lat, lon, resolution
                )
                if dggrid_id not in dggrid_ids_band_values:
                    vals = raster_data[:, int(row), int(col)]
                    dggrid_ids_band_values[dggrid_id] = [
                        (v.item() if hasattr(v, "item") else v) for v in vals
                    ]
            except Exception:
                continue

    # Build GeoDataFrame as the base
    properties = []
    for dggrid_id, band_values in tqdm(
        dggrid_ids_band_values.items(),
        desc="Converting raster to DGGRID",
        unit=" cells",
    ):
        try:
            # Convert zone to geometry using dggrid2geo
            cell_polygon = dggrid2geo(
                dggrid_instance,
                dggs_type,
                dggrid_id,
                resolution,
                split_antimeridian=split_antimeridian,
                aggregate=aggregate,
            )

            # Get cell metrics
            if isinstance(cell_polygon, gpd.GeoDataFrame) and not cell_polygon.empty:
                cell_geom = cell_polygon.iloc[0].geometry
                centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
                    geodesic_dggs_metrics(cell_geom, len(cell_geom.exterior.coords) - 1)
                )
            else:
                # Fallback metrics if geometry conversion fails
                centroid_lat, centroid_lon, avg_edge_len, cell_area, cell_perimeter = (
                    0,
                    0,
                    0,
                    0,
                    0,
                )

            base_props = {
                f"dggrid_{dggs_type}": dggrid_id,
                "resolution": resolution,
                "center_lat": centroid_lat,
                "center_lon": centroid_lon,
                "avg_edge_len": avg_edge_len,
                "cell_area": cell_area,
                "cell_perimeter": cell_perimeter,
                "geometry": cell_geom
                if isinstance(cell_polygon, gpd.GeoDataFrame) and not cell_polygon.empty
                else None,
            }
            band_properties = {
                f"band_{i + 1}": band_values[i] for i in range(band_count)
            }
            base_props.update(band_properties)
            properties.append(base_props)
        except Exception:
            continue

    # Build GeoDataFrame
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2dggrid" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to QTM Module

This module provides functionality to convert raster data to QTM (Quaternary Triangular Mesh) DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_qtm_resolution(raster_path)

Automatically determine the optimal QTM resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate QTM resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal QTM resolution level

Examples

cell_size, resolution = get_nearest_qtm_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2qtm.py 
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
def get_nearest_qtm_resolution(raster_path):
    """
    Automatically determine the optimal QTM resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate QTM resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal QTM resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_qtm_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = qtm_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2qtm(raster_path, resolution=None, output_format='gpd')

Convert raster data to QTM DGGS format.

Converts raster data to QTM (Quaternary Triangular Mesh) DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a QTM cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional QTM resolution level. If None, automatically determined based on raster pixel size. Valid range: 1-30. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a QTM cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2qtm("data.tif") print(f"Converted {len(result)} QTM cells")

Convert with specific resolution

result = raster2qtm("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2qtm("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2qtm.py 
 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
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
181
def raster2qtm(raster_path, resolution=None, output_format="gpd"):
    """
    Convert raster data to QTM DGGS format.

    Converts raster data to QTM (Quaternary Triangular Mesh) DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a QTM cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        QTM resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 1-30.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a QTM cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2qtm("data.tif")
    >>> print(f"Converted {len(result)} QTM cells")

    >>> # Convert with specific resolution
    >>> result = raster2qtm("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2qtm("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest qtm resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_qtm_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest QTM resolution determined: {resolution}")
    else:
        resolution = validate_qtm_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per QTM cell
    qtm_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            qtm_id = latlon2qtm(lat, lon, resolution)
            if qtm_id not in qtm_band_values:
                vals = raster_data[:, int(row), int(col)]
                qtm_band_values[qtm_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    # Always convert to GeoDataFrame for output
    # Build GeoDataFrame
    properties = []
    for qtm_id, band_values in tqdm(
        qtm_band_values.items(), desc="Converting raster to QTM", unit=" cells"
    ):
        cell_polygon = qtm2geo(qtm_id)
        base_props = {"qtm": qtm_id, "geometry": cell_polygon}
        band_props = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_props)
        properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2qtm" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

raster2qtm_cli()

Command line interface for raster2qtm conversion

Source code in vgrid/conversion/raster2dggs/raster2qtm.py 
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
def raster2qtm_cli():
    """Command line interface for raster2qtm conversion"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to QTM DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")
    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        default=None,
        help=f"QTM resolution [{min_res}..{max_res}]",
    )
    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )

    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format

    if not os.path.exists(raster):
        print(f"Error: The file {raster} does not exist.")
        return

    result = raster2qtm(raster, resolution, output_format)
    if output_format in STRUCTURED_FORMATS:
        print(result)

Raster to OLC Module

This module provides functionality to convert raster data to OLC (Open Location Code) DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_olc_resolution(raster_path)

Automatically determine the optimal OLC resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate OLC resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal OLC resolution level

Examples

cell_size, resolution = get_nearest_olc_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 8

Source code in vgrid/conversion/raster2dggs/raster2olc.py 
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
def get_nearest_olc_resolution(raster_path):
    """
    Automatically determine the optimal OLC resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate OLC resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal OLC resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_olc_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 8
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    # Find the nearest s2 resolution by comparing the pixel size to the s2 edge lengths
    nearest_resolution = None
    min_diff = float("inf")

    for res in olc_res:
        _, _, avg_area, _ = olc_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2olc(raster_path, resolution=None, output_format='gpd')

Convert raster data to OLC DGGS format.

Converts raster data to OLC (Open Location Code) DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to an OLC cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional OLC resolution level. If None, automatically determined based on raster pixel size. Valid values: [2, 4, 6, 8, 10, 11, 12, 13, 14, 15]. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents an OLC cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2olc("data.tif") print(f"Converted {len(result)} OLC cells")

Convert with specific resolution

result = raster2olc("data.tif", resolution=8)

Convert to GeoJSON file

result = raster2olc("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2olc.py 
 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
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
181
182
183
184
185
def raster2olc(raster_path, resolution=None, output_format="gpd"):
    """
    Convert raster data to OLC DGGS format.

    Converts raster data to OLC (Open Location Code) DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to an OLC cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        OLC resolution level. If None, automatically determined based on raster pixel size.
        Valid values: [2, 4, 6, 8, 10, 11, 12, 13, 14, 15].
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents an OLC cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2olc("data.tif")
    >>> print(f"Converted {len(result)} OLC cells")

    >>> # Convert with specific resolution
    >>> result = raster2olc("data.tif", resolution=8)

    >>> # Convert to GeoJSON file
    >>> result = raster2olc("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest olc resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_olc_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest OLC resolution determined: {resolution}")
    else:
        resolution = validate_olc_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per OLC cell
    olc_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            olc_id = latlon2olc(lat, lon, resolution)
            if olc_id not in olc_band_values:
                vals = raster_data[:, int(row), int(col)]
                olc_band_values[olc_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for olc_id, band_values in tqdm(
        olc_band_values.items(), desc="Converting raster to OLC", unit=" cells"
    ):
        cell_polygon = olc2geo(olc_id)
        min_lon, min_lat, max_lon, max_lat = cell_polygon.bounds
        cell_polygon = Polygon(
            [
                [min_lon, min_lat],
                [max_lon, min_lat],
                [max_lon, max_lat],
                [min_lon, max_lat],
                [min_lon, min_lat],
            ]
        )
        base_props = {"olc": olc_id, "geometry": cell_polygon}
        band_props = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_props)
        properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2olc" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

raster2olc_cli()

Command line interface for raster to OLC conversion

Source code in vgrid/conversion/raster2dggs/raster2olc.py 
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def raster2olc_cli():
    """Command line interface for raster to OLC conversion"""
    parser = argparse.ArgumentParser(
        description="Convert Raster in Geographic CRS to OLC/ Google Plus Code DGGS"
    )
    parser.add_argument("-raster", type=str, required=True, help="Raster file path")

    parser.add_argument(
        "-r",
        "--resolution",
        type=int,
        required=False,
        choices=olc_res,
        default=None,
        help="OLC resolution",
    )

    parser.add_argument(
        "-f",
        "--output_format",
        type=str,
        choices=OUTPUT_FORMATS,
        default="gpd",
    )

    args = parser.parse_args()
    raster = args.raster
    resolution = args.resolution
    output_format = args.output_format

    if not os.path.exists(raster):
        print(f"Error: The file {raster} does not exist.")
        return

    result = raster2olc(raster, resolution, output_format)
    if output_format in STRUCTURED_FORMATS:
        print(result)

Raster to Geohash Module

This module provides functionality to convert raster data to Geohash DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_geohash_resolution(raster_path)

Automatically determine the optimal Geohash resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate Geohash resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal Geohash resolution level

Examples

cell_size, resolution = get_nearest_geohash_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2geohash.py 
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
def get_nearest_geohash_resolution(raster_path):
    """
    Automatically determine the optimal Geohash resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate Geohash resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal Geohash resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_geohash_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = geohash_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res
    return cell_size, nearest_resolution

raster2geohash(raster_path, resolution=None, output_format='gpd')

Convert raster data to Geohash DGGS format.

Converts raster data to Geohash DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a Geohash cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional Geohash resolution level. If None, automatically determined based on raster pixel size. Valid range: 1-12. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a Geohash cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2geohash("data.tif") print(f"Converted {len(result)} Geohash cells")

Convert with specific resolution

result = raster2geohash("data.tif", resolution=5)

Convert to GeoJSON file

result = raster2geohash("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2geohash.py 
 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
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
def raster2geohash(raster_path, resolution=None, output_format="gpd"):
    """
    Convert raster data to Geohash DGGS format.

    Converts raster data to Geohash DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a Geohash cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        Geohash resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 1-12.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a Geohash cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2geohash("data.tif")
    >>> print(f"Converted {len(result)} Geohash cells")

    >>> # Convert with specific resolution
    >>> result = raster2geohash("data.tif", resolution=5)

    >>> # Convert to GeoJSON file
    >>> result = raster2geohash("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest geohash resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_geohash_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest Geohash resolution determined: {resolution}")
    else:
        resolution = validate_geohash_resolution(resolution)

    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per Geohash cell
    geohash_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            geohash_id = latlon2geohash(lat, lon, resolution)
            if geohash_id not in geohash_band_values:
                vals = raster_data[:, int(row), int(col)]
                geohash_band_values[geohash_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for geohash_id, band_values in tqdm(
        geohash_band_values.items(), desc="Converting raster to Geohash", unit=" cells"
    ):
        cell_polygon = geohash2geo(geohash_id)
        base_props = {"geohash": geohash_id, "geometry": cell_polygon}
        band_props = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_props)
        properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2geohash" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to Tilecode Module

This module provides functionality to convert raster data to Tilecode DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_tilecode_resolution(raster_path)

Automatically determine the optimal Tilecode resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate Tilecode resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal Tilecode resolution level

Examples

cell_size, resolution = get_nearest_tilecode_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2tilecode.py 
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
def get_nearest_tilecode_resolution(raster_path):
    """
    Automatically determine the optimal Tilecode resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate Tilecode resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal Tilecode resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_tilecode_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = tilecode_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2tilecode(raster_path, resolution=None, output_format='gpd')

Convert raster data to Tilecode DGGS format.

Converts raster data to Tilecode DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a Tilecode cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional Tilecode resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-26. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a Tilecode cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2tilecode("data.tif") print(f"Converted {len(result)} Tilecode cells")

Convert with specific resolution

result = raster2tilecode("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2tilecode("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2tilecode.py 
 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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def raster2tilecode(raster_path, resolution=None, output_format="gpd"):
    """
    Convert raster data to Tilecode DGGS format.

    Converts raster data to Tilecode DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a Tilecode cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        Tilecode resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-26.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a Tilecode cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2tilecode("data.tif")
    >>> print(f"Converted {len(result)} Tilecode cells")

    >>> # Convert with specific resolution
    >>> result = raster2tilecode("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2tilecode("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest tilecode resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_tilecode_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest Tilecode resolution determined: {resolution}")
    else:
        resolution = validate_tilecode_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per Tilecode cell
    tilecode_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            tilecode_id = tilecode.latlon2tilecode(lat, lon, resolution)
            if tilecode_id not in tilecode_band_values:
                vals = raster_data[:, int(row), int(col)]
                tilecode_band_values[tilecode_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for tilecode_id, band_values in tqdm(
        tilecode_band_values.items(),
        desc="Converting raster to Tilecode",
        unit=" cells",
    ):
        match = re.match(r"z(\d+)x(\d+)y(\d+)", tilecode_id)
        if match:
            z = int(match.group(1))
            x = int(match.group(2))
            y = int(match.group(3))
            bounds = mercantile.bounds(x, y, z)
            if bounds:
                min_lat, min_lon = bounds.south, bounds.west
                max_lat, max_lon = bounds.north, bounds.east
                cell_polygon = Polygon(
                    [
                        [min_lon, min_lat],
                        [max_lon, min_lat],
                        [max_lon, max_lat],
                        [min_lon, max_lat],
                        [min_lon, min_lat],
                    ]
                )
                base_props = {"tilecode": tilecode_id, "geometry": cell_polygon}
                band_props = {
                    f"band_{i + 1}": band_values[i] for i in range(band_count)
                }
                base_props.update(band_props)
                properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2tilecode" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)

Raster to Quadkey Module

This module provides functionality to convert raster data to Quadkey DGGS format with automatic resolution determination and multi-band support.

Key Functions

get_nearest_quadkey_resolution(raster_path)

Automatically determine the optimal Quadkey resolution for a given raster.

Analyzes the raster's pixel size and determines the most appropriate Quadkey resolution that best matches the raster's spatial resolution.

Parameters

raster_path : str Path to the raster file to analyze.

Returns

tuple A tuple containing (cell_size, resolution) where: - cell_size: The calculated cell size in square meters - resolution: The optimal Quadkey resolution level

Examples

cell_size, resolution = get_nearest_quadkey_resolution("data.tif") print(f"Cell size: {cell_size} m², Resolution: {resolution}") Cell size: 1000000.0 m², Resolution: 5

Source code in vgrid/conversion/raster2dggs/raster2quadkey.py 
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
def get_nearest_quadkey_resolution(raster_path):
    """
    Automatically determine the optimal Quadkey resolution for a given raster.

    Analyzes the raster's pixel size and determines the most appropriate Quadkey resolution
    that best matches the raster's spatial resolution.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to analyze.

    Returns
    -------
    tuple
        A tuple containing (cell_size, resolution) where:
        - cell_size: The calculated cell size in square meters
        - resolution: The optimal Quadkey resolution level

    Examples
    --------
    >>> cell_size, resolution = get_nearest_quadkey_resolution("data.tif")
    >>> print(f"Cell size: {cell_size} m², Resolution: {resolution}")
    Cell size: 1000000.0 m², Resolution: 5
    """
    with rasterio.open(raster_path) as src:
        transform = src.transform
        crs = src.crs
        pixel_width = transform.a
        pixel_height = -transform.e
        cell_size = pixel_width * pixel_height

        if crs.is_geographic:
            # Latitude of the raster center
            center_latitude = (src.bounds.top + src.bounds.bottom) / 2
            # Convert degrees to meters
            meter_per_degree_lat = 111_320  # Roughly 1 degree latitude in meters
            meter_per_degree_lon = meter_per_degree_lat * cos(radians(center_latitude))

            pixel_width_m = pixel_width * meter_per_degree_lon
            pixel_height_m = pixel_height * meter_per_degree_lat
            cell_size = pixel_width_m * pixel_height_m

    min_diff = float("inf")
    nearest_resolution = min_res

    for res in range(min_res, max_res + 1):
        _, _, avg_area, _ = quadkey_metrics(res)
        if avg_area < MIN_CELL_AREA:
            break
        diff = abs(avg_area - cell_size)
        # If the difference is smaller than the current minimum, update the nearest resolution
        if diff < min_diff:
            min_diff = diff
            nearest_resolution = res

    return cell_size, nearest_resolution

raster2quadkey(raster_path, resolution=None, output_format='gpd')

Convert raster data to Quadkey DGGS format.

Converts raster data to Quadkey DGGS format with automatic resolution determination and multi-band support. Each pixel is assigned to a Quadkey cell and the first sample value per cell is preserved.

Parameters

raster_path : str Path to the raster file to convert. resolution : int, optional Quadkey resolution level. If None, automatically determined based on raster pixel size. Valid range: 0-29. output_format : str, default "gpd" Output format. Options: - "gpd": Returns GeoPandas GeoDataFrame (default) - "csv": Returns CSV file path - "geojson": Returns GeoJSON file path - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict - "parquet": Returns Parquet file path - "shapefile"/"shp": Returns Shapefile file path - "gpkg"/"geopackage": Returns GeoPackage file path

Returns

geopandas.GeoDataFrame or str or dict The converted data in the specified format. Each row represents a Quadkey cell with geometry and band values from the original raster.

Examples

Convert with automatic resolution

result = raster2quadkey("data.tif") print(f"Converted {len(result)} Quadkey cells")

Convert with specific resolution

result = raster2quadkey("data.tif", resolution=10)

Convert to GeoJSON file

result = raster2quadkey("data.tif", output_format="geojson") print(f"Saved to: {result}")

Source code in vgrid/conversion/raster2dggs/raster2quadkey.py 
 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
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
181
182
183
184
185
186
187
188
189
190
191
192
193
194
def raster2quadkey(raster_path, resolution=None, output_format="gpd"):
    """
    Convert raster data to Quadkey DGGS format.

    Converts raster data to Quadkey DGGS format with automatic resolution
    determination and multi-band support. Each pixel is assigned to a Quadkey cell and
    the first sample value per cell is preserved.

    Parameters
    ----------
    raster_path : str
        Path to the raster file to convert.
    resolution : int, optional
        Quadkey resolution level. If None, automatically determined based on raster pixel size.
        Valid range: 0-29.
    output_format : str, default "gpd"
        Output format. Options:
        - "gpd": Returns GeoPandas GeoDataFrame (default)
        - "csv": Returns CSV file path
        - "geojson": Returns GeoJSON file path
        - "geojson_dict": Returns GeoJSON FeatureCollection as Python dict
        - "parquet": Returns Parquet file path
        - "shapefile"/"shp": Returns Shapefile file path
        - "gpkg"/"geopackage": Returns GeoPackage file path

    Returns
    -------
    geopandas.GeoDataFrame or str or dict
        The converted data in the specified format. Each row represents a Quadkey cell
        with geometry and band values from the original raster.

    Examples
    --------
    >>> # Convert with automatic resolution
    >>> result = raster2quadkey("data.tif")
    >>> print(f"Converted {len(result)} Quadkey cells")

    >>> # Convert with specific resolution
    >>> result = raster2quadkey("data.tif", resolution=10)

    >>> # Convert to GeoJSON file
    >>> result = raster2quadkey("data.tif", output_format="geojson")
    >>> print(f"Saved to: {result}")
    """
    # Step 1: Determine the nearest quadkey resolution if none is provided
    if resolution is None:
        cell_size, resolution = get_nearest_quadkey_resolution(raster_path)
        print(f"Cell size: {cell_size} m2")
        print(f"Nearest Quadkey resolution determined: {resolution}")
    else:
        resolution = validate_quadkey_resolution(resolution)
    # Open the raster file to get metadata and data
    with rasterio.open(raster_path) as src:
        raster_data = src.read()  # Read all bands
        transform = src.transform
        width, height = src.width, src.height
        band_count = src.count  # Number of bands in the raster

    # Collect band values during the pixel scan, storing the first sample per Quadkey cell
    quadkey_band_values = {}
    for row in range(height):
        for col in range(width):
            lon, lat = transform * (col, row)
            quadkey_id = tilecode.latlon2quadkey(lat, lon, resolution)
            if quadkey_id not in quadkey_band_values:
                vals = raster_data[:, int(row), int(col)]
                quadkey_band_values[quadkey_id] = [
                    (v.item() if hasattr(v, "item") else v) for v in vals
                ]

    properties = []
    for quadkey_id, band_values in tqdm(
        quadkey_band_values.items(), desc="Converting raster to Quadkey", unit=" cells"
    ):
        tile = mercantile.quadkey_to_tile(quadkey_id)
        z = tile.z
        x = tile.x
        y = tile.y
        bounds = mercantile.bounds(x, y, z)
        min_lat, min_lon = bounds.south, bounds.west
        max_lat, max_lon = bounds.north, bounds.east
        cell_polygon = Polygon(
            [
                [min_lon, min_lat],
                [max_lon, min_lat],
                [max_lon, max_lat],
                [min_lon, max_lat],
                [min_lon, min_lat],
            ]
        )
        base_props = {"quadkey": quadkey_id, "geometry": cell_polygon}
        band_props = {f"band_{i + 1}": band_values[i] for i in range(band_count)}
        base_props.update(band_props)
        properties.append(base_props)
    gdf = gpd.GeoDataFrame(properties, geometry="geometry", crs="EPSG:4326")

    # Use centralized output utility
    base_name = os.path.splitext(os.path.basename(raster_path))[0]
    output_name = f"{base_name}2quadkey" if output_format is not None else None
    return convert_to_output_format(gdf, output_format, output_name)