计算矢量中多边形总面积、栅格范围裁剪为小正方形

目录

计算矢量中多边形总面积

读取栅格，将栅格裁剪为256*256个正方形小栅格，舍弃边缘。将小栅格范围转矢量，输出矢量

---

计算矢量中多边形总面积

from osgeo import ogr, osr

def calculate_area(file_path):
    # 打开矢量文件
    dataSource = ogr.Open(file_path)
    if dataSource is None:
        print("Could not open file")
        return

    layer = dataSource.GetLayer()

    # 获取层的空间参考系
    source_srs = layer.GetSpatialRef()

    # 创建目标空间参考系（WGS 84 / World Mercator, EPSG:3395）
    target_srs = osr.SpatialReference()
    target_srs.ImportFromEPSG(3395)  # EPSG:3395 是米单位的投影坐标系

    # 创建坐标变换
    transform = osr.CoordinateTransformation(source_srs, target_srs)

    total_area = 0.0

    # 遍历每个要素并计算面积
    for feature in layer:
        geom = feature.GetGeometryRef()
        geom.Transform(transform)  # 转换到目标空间参考系
        area = geom.GetArea()  # 获取面积，单位为平方米
        total_area += area

    # 将面积从平方米转换为平方千米
    total_area_km2 = total_area / 1e6

    print(f"Total Area: {total_area_km2:.2f} square kilometers")

# 替换为你的矢量文件路径
file_path = 'unet/2.shp'
calculate_area(file_path)

读取栅格，将栅格裁剪为256*256个正方形小栅格，舍弃边缘。将小栅格范围转矢量，输出矢量

import os
from osgeo import gdal, ogr, osr

def raster_to_polygons(input_raster, output_shapefile, tile_size):
    # Open the input raster file
    raster = gdal.Open(input_raster)
    if not raster:
        raise FileNotFoundError(f"Unable to open raster file {input_raster}")
    
    # Get raster georeference info
    transform = raster.GetGeoTransform()
    pixel_width = transform[1]
    pixel_height = transform[5]
    origin_x = transform[0]
    origin_y = transform[3]
    
    # Get raster size
    cols = raster.RasterXSize
    rows = raster.RasterYSize
    
    # Calculate the number of tiles in x and y directions
    num_tiles_x = cols // tile_size
    num_tiles_y = rows // tile_size
    
    # Define the spatial reference
    srs = osr.SpatialReference()
    srs.ImportFromWkt(raster.GetProjection())
    
    # Create the output shapefile
    driver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(output_shapefile):
        driver.DeleteDataSource(output_shapefile)
    out_ds = driver.CreateDataSource(output_shapefile)
    out_layer = out_ds.CreateLayer("tiles", srs, geom_type=ogr.wkbPolygon)
    
    # Add fields to the shapefile
    id_field = ogr.FieldDefn("ID", ogr.OFTInteger)
    out_layer.CreateField(id_field)
    
    tile_id = 0
    
    # Loop through the tiles and create polygons
    for i in range(num_tiles_x):
        for j in range(num_tiles_y):
            min_x = origin_x + i * tile_size * pixel_width
            max_x = min_x + tile_size * pixel_width
            min_y = origin_y + j * tile_size * pixel_height
            max_y = min_y + tile_size * pixel_height
            
            # Create polygon geometry
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(min_x, min_y)
            ring.AddPoint(max_x, min_y)
            ring.AddPoint(max_x, max_y)
            ring.AddPoint(min_x, max_y)
            ring.AddPoint(min_x, min_y)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)
            
            # Create feature and set geometry and attributes
            feature = ogr.Feature(out_layer.GetLayerDefn())
            feature.SetGeometry(poly)
            feature.SetField("ID", tile_id)
            out_layer.CreateFeature(feature)
            feature = None
            
            tile_id += 1
    
    # Close the shapefile
    out_ds = None
    
    print(f"{tile_id} tiles created and saved to {output_shapefile}")

# Example usage
input_raster = '指定范围真实标签.tif'
output_shapefile = '1.shp'
tile_size = 256  # Size of each tile in pixels

raster_to_polygons(input_raster, output_shapefile, tile_size)

要素转栅格