tt.py
1.1 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
# coding=utf-8
#author: 4N
#createtime: 2021/7/15
#email: nheweijun@sina.com
from osgeo import gdal
from osgeo import gdal, gdalconst
from osgeo import ogr
rasterFile = 'F:/**0416.dat' # 原影像
shpFile = 'F:/**小麦.shp' # 裁剪矩形
dataset = gdal.Open(rasterFile, gdalconst.GA_ReadOnly)
geo_transform = dataset.GetGeoTransform()
cols = dataset.RasterXSize # 列数
rows = dataset.RasterYSize # 行数
x_min = geo_transform[0]
y_min = geo_transform[3]
pixel_width = geo_transform[1]
shp = ogr.Open(shpFile, 0)
m_layer = shp.GetLayerByIndex(0)
target_ds = gdal.GetDriverByName('MEM').Create("", xsize=cols, ysize=rows, bands=1,
eType=gdal.GDT_Byte)
target_ds.SetGeoTransform(geo_transform)
target_ds.SetProjection(dataset.GetProjection())
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(0)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], m_layer) # 跟shp字段给栅格像元赋值
# gdal.RasterizeLayer(target_ds, [1], m_layer) # 多边形内像元值的全是255
del dataset
del target_ds
shp.Release()