Растеризация полигона с использованием координат в градусах в GDAL Python - PullRequest
0 голосов
/ 15 февраля 2020

Я пытаюсь растеризовать полигон из координат в градусах, используя GDAL и Python. Пока что он работает с использованием координат в метрах, но они зависят от проекции системы координат. Я мог бы преобразовать координаты в файле формы из градусов в метры до растеризации многоугольника, что решило бы проблему, но есть ли способ вызвать RasterizeLayer, используя координаты градусов?

def create_polygon_shapefile(coords, shp_path):
    # create polygon object
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for coord in coords:
        ring.AddPoint(coord[0], coord[1])
    polygon = ogr.Geometry(ogr.wkbPolygon)
    polygon.AddGeometry(ring)
    driver = ogr.GetDriverByName('ESRI Shapefile')
    ds = driver.CreateDataSource(shp_path)
    # create layer and add polygon
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)
    layer = ds.CreateLayer('', srs, ogr.wkbPolygon)
    fieldDefn_ = ogr.FieldDefn('id', ogr.OFTInteger)
    layer.CreateField(fieldDefn_)
    featureDefn = layer.GetLayerDefn()
    feature = ogr.Feature(featureDefn)
    feature.SetGeometry(polygon)
    feature.SetField('id', 1)
    layer.CreateFeature(feature)


def rasterize_polygon(shp_file, target_file, w, h, gt):
    source_ds = ogr.Open(shp_file)
    layer = source_ds.GetLayer()

    # create target raster
    driver = gdal.GetDriverByName('GTiff')
    target_ds = driver.Create(target_file, w, h, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform(gt)

    # rasterize polygon onto target raster
    gdal.RasterizeLayer(target_ds, [1], layer)
    target_ds.FlushCache()
    return target_ds.GetRasterBand(1).ReadAsArray()

# a small example:
w, h = 20, 20
gt = (0, 1, 0, 0, 0, 1)
coords = [(10, 2), (18, 14), (3, 16), (8, 12)]
create_polygon_shapefile(coords, 'temp.shp')

for row in rasterize_polygon('temp.shp', 'target.tif', w, h, gt):
    for e in row:
        if e == 0:
            print(' ', end='')
        else:
            print('#', end='')
    print()
...