GDAL学习记录

一、osgeo常识

1. gdal: 处理栅格

   • 重采样方法：gdal.GRA_*
   • 数据类型：gdal.GDT_*
   • 关键类：gdal.Dataset, gdal.Driver, gdal.Band

2. osr:【proj = osr.SpatialReference().SetWellKnownGeogCS(‘WGS84’)】

   • wkt: ‘GEOGCS[“WGS 84”,…]’
   • proj4: ‘+proj=longlat +datum=WGS84 +no_defs’
   • epsg: ‘4326’,【proj.GetAuthorityCode(None)】
     • 投影坐标系：32601-32660，Central Meridian = (code-32601)*6-117

3. ogr: 处理矢量

   • 几何类型: ogr.wkb*

   • 字段值类型：ogr.OFT*

     1	2	3	4	5	6	7
     ogr.DataSource	GetDriver	GetLayer
     GetLayerCount	ogr.Driver	ogr.Layer	GetFeature	GetLayerDefn	GetGeomType	GetSpatialRef
     		GetExtent GetFeatureCount	ogr.Feature	ogr.FeatureDefn	ogr.wkb*	osr.SpatialReference
     ogr.Feature	GetFieldDefn	GetFieldType	GetGeometryRef	–	ogr.FeatureDefn	GetFieldDefn
     GetFieldCount GetField	ogr.FieldDefn	ogr.OFT	ogr.Geometry	GetSpatialReference	GetFieldCount	ogr.FieldDefn
     –	–	–	GetGeometryName GetGeometryCount GetPointCount GetEnvelope	osr.SpatialReference		GetName GetWidth GetType GetPrecision

二、Warp和WarpOptions

gdal.WarpOptions(): Create a WarpOptions() object that can be passed to gdal.Warp()

• format — 输出格式 (“GTiff”, etc…)

• outputType — output type (gdalconst.GDT_Byte, etc…)

• srcNodata, dstNodata — source(output) nodata value(s)

• multithread — whether to multithread computation and I/O operations

• creationOptions — list of creation options, eg: [‘COMPRESS=LZW’,“TILED=True”]

  • TILED=YES/NO：指定输出数据集是否为分块格式。
  • COMPRESS=[ “JPEG”, “LZW”, “PACKBITS”, “DEFLATE”, “ZSTD”]：指定输出数据集的压缩格式；
  • BLOCKXSIZE=SIZE 和 BLOCKYSIZE=SIZE：指定输出数据集的块大小。
  • NUM_THREADS=[<integer>|NUM_CPUS]:

• warpOptions — list of warping options, eg: [‘NUM_THREADS=ALL_CPUS’]

  • WRITE_FLUSH=YES/NO: 强制在处理每个块之后刷新到数据磁盘
  • NUM_THREADS=integeror ALL_CPUS: 用于并行化warp计算的线程数

• coordinateOperation – coordinate operation as a PROJ string or WKT string

• copyMetadata — whether to copy source metadata

• overviewLevel — To specify which overview level of source files must be used

• 影像裁剪

  • cutlineDSName — 分割线数据名
  • cutlineLayer — 分割线图层名
  • cutlineWhere — 分割线where字句
  • cutlineSQL — 分割线SQL声明
  • cutlineBlend — 以像素为单位的分割线混合距离
  • cropToCutline — 是否使用分割线范围作为输出界限

• 影像镶嵌

• 影像重采样: 【xRes, yRes, srcSRS, dstSRS, resampleAlg】

• 影像重投影：【srcSRS, dstSRS, resampleAlg】

• 几何校正/正射校正

  • rpc — whether to use RPC transformer
  • polynomialOrder — order of polynomial GCP interpolation
  • transformerOptions — list of transformer options, eg:[‘RPC_DEM=’ + dem]

1、影像裁剪

src = gdal.Open(in_tif)
warpOptions = gdal.WarpOptions(format='GTiff', cutlineDSName=in_shp_r, cropToCutline=True)
gdal.Warp(out_tif_clip_r, src, options=warpOptions)
gdal.Warp(out_tif_clip_nr, src, format='GTiff', cutlineDSName=in_shp_nr, cropToCutline=True)

2、影像重投影、重采样

src_clip = gdal.Open(out_tif_clip_r)
gdal.Warp(out_tif_reproject, src_clip, dstSRS='EPSG:32650', resampleAlg=gdal.GRA_Bilinear)

out_tif_reproject = gdal.Open(out_tif_reproject)
gdal.Warp(out_tif_resample, out_tif_reproject, format='GTiff', xRes=10, yRes=10, resampleAlg=gdal.GRA_Bilinear)

3、影像镶嵌

s2a1, s2a2 = r"F:\RSdataLocal\S2A_20231024\T50SMJ.tif", r"F:\RSdataLocal\S2A_20231024\T50TMK.tif"
out_tif_mosaic = r"F:\RSdataLocal\S2A_20231024\T50_mosaic.tif"
in_proj = gdal.Open(s2a1).GetProjection()
gdal.Warp(out_tif_mosaic, [s2a1, s2a2], format='GTiff', srcSRS=in_proj, copyMetadata=True)

fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 5))
show_rgb(out_tif_clip_r, ax1, vmax=3000), show_rgb(out_tif_clip_nr, ax2, vmax=3000)