Geometry2Raster.py
1.8 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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
# coding=utf-8
#author: 4N
#createtime: 2021/9/29
#email: nheweijun@sina.com
from osgeo.ogr import Layer,Geometry
from osgeo import gdal,ogr
class Geometry2Raster:
def __init__(self):
pass
@classmethod
def convert(cls,geometry:Geometry,bbox, xy_res):
'''
根据Geometry和bbox将空间对象栅格化
:param geometry: 空间对象
:param bbox: 栅格化范围
:param xy_res: 分辨率
:return:
'''
#计算像素大小
x_min,y_min,x_max, y_max = bbox
pixel_size = (x_max - x_min) / xy_res
#创建内存栅格对象
target_ds = gdal.GetDriverByName('MEM').Create('', xy_res, xy_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
#创建内存图层对象，因为gdal只提供栅格化图层的方法
layer_ds = ogr.GetDriverByName('Memory').CreateDataSource('')
mem_layer= layer_ds.CreateLayer("mem", geom_type=ogr.wkbUnknown)
feature_defn = mem_layer.GetLayerDefn()
feature = ogr.Feature(feature_defn)
feature.SetGeometry(geometry)
mem_layer.CreateFeature(feature)
#栅格化
gdal.RasterizeLayer(target_ds, [1], mem_layer, burn_values=[1])
data = band.ReadAsArray()
del target_ds
del layer_ds
return data
if __name__ == '__main__':
# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(111.62, 29)
ring.AddPoint(111.62, 29.03)
ring.AddPoint(111.727, 29.03)
ring.AddPoint(111.727, 29)
ring.AddPoint(111.62, 29)
# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)
poly.ExportToWkt()