Я пытаюсь растеризовать полигон из координат в градусах, используя 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()