from raster_conversion import RasterConverter
import os

input_file = "path/to/input_raster.tif"
output_file = os.path.join(os.path.dirname(input_file), "output_boundary.shp")

if __name__ == '__main__':
    converter = RasterConverter()
    converter.process(input_file, output_file)

Core Implementation

Class Definition

import os
import time
import gc
from multiprocessing import Pool, cpu_count
from osgeo import gdal, ogr, osr
import numpy as np

class RasterConverter:
    def __init__(self, block_size=512, memory_limit=256, use_parallel=True):
        gdal.UseExceptions()
        self.block_size = block_size
        self.memory_limit = memory_limit
        self.use_parallel = use_parallel
        self.core_count = cpu_count()

Raster Information Extraction

    def get_raster_metadata(self, raster_path):
        dataset = gdal.Open(raster_path, gdal.GA_ReadOnly)
        if not dataset:
            raise ValueError(f"Cannot open raster file: {raster_path}")
        
        metadata = {
            'dimensions': (dataset.RasterXSize, dataset.RasterYSize),
            'band_count': dataset.RasterCount,
            'transform': dataset.GetGeoTransform(),
            'projection': dataset.GetProjection(),
            'data_type': dataset.GetRasterBand(1).DataType,
            'no_data': dataset.GetRasterBand(1).GetNoDataValue()
        }
        
        pixel_count = metadata['dimensions'][0] * metadata['dimensions'][1]
        bytes_per_pixel = gdal.GetDataTypeSize(metadata['data_type']) // 8
        metadata['estimated_memory'] = (pixel_count * bytes_per_pixel) / (1024 * 1024)
        
        dataset = None
        return metadata

Processing Strateyg Selection

    def needs_block_processing(self, metadata):
        print(f"Estimated memory requirement: {metadata['estimated_memory']:.1f} MB")
        return (metadata['estimated_memory'] > self.memory_limit or
                metadata['dimensions'][0] > self.block_size * 4 or
                metadata['dimensions'][1] > self.block_size * 4)

Mask Generation

    def generate_mask(self, raster_path, no_data=None):
        dataset = gdal.Open(raster_path, gdal.GA_ReadOnly)
        if not dataset:
            raise ValueError(f"Cannot open raster file: {raster_path}")
        
        width, height = dataset.RasterXSize, dataset.RasterYSize
        transform = dataset.GetGeoTransform()
        projection = dataset.GetProjection()
        
        band = dataset.GetRasterBand(1)
        if no_data is None:
            no_data = band.GetNoDataValue()
        
        data = band.ReadAsArray()
        
        if no_data is not None:
            mask = np.where(np.isclose(data, no_data, equal_nan=True), 0, 1)
        else:
            mask = np.where(np.isnan(data), 0, 1)
        
        mask = mask.astype(np.uint8)
        dataset = None
        
        return mask, transform, projection, width, height

Block Processing Implementation

    def generate_mask_blocks(self, raster_path, no_data=None):
        print("Using block processing mode...")
        start = time.time()
        
        dataset = gdal.Open(raster_path, gdal.GA_ReadOnly)
        width, height = dataset.RasterXSize, dataset.RasterYSize
        transform = dataset.GetGeoTransform()
        projection = dataset.GetProjection()
        
        band = dataset.GetRasterBand(1)
        if no_data is None:
            no_data = band.GetNoDataValue()
        
        mask = np.zeros((height, width), dtype=np.uint8)
        
        x_blocks = (width + self.block_size - 1) // self.block_size
        y_blocks = (height + self.block_size - 1) // self.block_size
        
        for y_block in range(y_blocks):
            for x_block in range(x_blocks):
                x_start = x_block * self.block_size
                y_start = y_block * self.block_size
                x_size = min(self.block_size, width - x_start)
                y_size = min(self.block_size, height - y_start)
                
                data_block = band.ReadAsArray(x_start, y_start, x_size, y_size)
                
                if no_data is not None:
                    mask_block = np.where(np.isclose(data_block, no_data, equal_nan=True), 0, 1)
                else:
                    mask_block = np.where(np.isnan(data_block), 0, 1)
                
                mask[y_start:y_start+y_size, x_start:x_start+x_size] = mask_block.astype(np.uint8)
                
                del data_block, mask_block
                gc.collect()
        
        dataset = None
        print(f"Block processing completed in {time.time() - start:.1f}s")
        return mask, transform, projection, width, height

Vetcorization Process

    def vectorize_mask(self, mask, transform, projection, output_path, format='ESRI Shapefile'):
        mem_driver = gdal.GetDriverByName('MEM')
        mem_dataset = mem_driver.Create('', mask.shape[1], mask.shape[0], 1, gdal.GDT_Byte)
        mem_dataset.SetGeoTransform(transform)
        mem_dataset.SetProjection(projection)
        
        mem_band = mem_dataset.GetRasterBand(1)
        mem_band.WriteArray(mask)
        mem_band.SetNoDataValue(0)
        
        if format == 'GeoJSON':
            driver = 'GeoJSON'
            if not output_path.endswith('.geojson'):
                output_path += '.geojson'
        else:
            driver = 'ESRI Shapefile'
            if not output_path.endswith('.shp'):
                output_path += '.shp'
        
        if os.path.exists(output_path):
            ogr.GetDriverByName(driver).DeleteDataSource(output_path)
        
        vector_ds = ogr.GetDriverByName(driver).CreateDataSource(output_path)
        srs = osr.SpatialReference()
        srs.ImportFromWkt(projection)
        
        layer = vector_ds.CreateLayer('boundary', srs, ogr.wkbPolygon)
        layer.CreateField(ogr.FieldDefn('VALUE', ogr.OFTInteger))
        
        gdal.Polygonize(mem_band, mem_band, layer, 0, [])
        
        mem_dataset = None
        vector_ds = None
        
        print(f"Vectorization complete: {output_path}")

Main Processing Method

    def process(self, input_raster, output_vector, format='ESRI Shapefile', no_data=None):
        try:
            print(f"Processing raster: {input_raster}")
            start = time.time()
            
            metadata = self.get_raster_metadata(input_raster)
            use_blocks = self.needs_block_processing(metadata)
            
            if use_blocks:
                mask, transform, projection, width, height = self.generate_mask_blocks(input_raster, no_data)
            else:
                mask, transform, projection, width, height = self.generate_mask(input_raster, no_data)
            
            temp_output = output_vector + '_temp'
            self.vectorize_mask(mask, transform, projection, temp_output, format)
            
            if format == 'ESRI Shapefile':
                temp_output += '.shp'
                for ext in ['.shp', '.shx', '.dbf', '.prj']:
                    if os.path.exists(temp_output.replace('.shp', ext)):
                        os.remove(temp_output.replace('.shp', ext))
            else:
                temp_output += '.geojson'
                if os.path.exists(temp_output):
                    os.remove(temp_output)
            
            del mask
            gc.collect()
            
            print(f"Processing completed in {time.time() - start:.1f}s")
        
        except Exception as e:
            print(f"Error during processing: {str(e)}")
            raise