一、测试gdal是否成功运行

import gdal
import os 
from matplotlib import pyplot as plt
import rasterio as ras
import numpy as np
from osgeo import gdal
import os
import pathlib
import rasterio
from rasterio.plot import show
import pandas as pd
import multiprocessing as mp
from itertools import repeat
#查看gdal版本
print(gdal.VersionInfo())
#打开栅格文件

path=r'F:\重庆数据\潼南-中分影像\32645490313\4903132021010100000000332239150100002.tif'
ds = gdal.Open(path)
if ds is None:
    print("打开失败！")
else:
#读取行列数
    print("读取行列数:")
    cols=ds.RasterXSize
    rows=ds.RasterYSize
    print("行数为：{} 列数为：{}".format(cols,rows))
#读取坐标系统
    print("读取坐标系统:")
    proj=ds.GetProjection()
    del ds
print(proj)
i=490313

2. 重投影和重采样

# %%time
"""

重投影时，重采样

重投影，投影过程中可以设置以下信息
gdal.Warp(
    xRes,yRes:两个方向上的分辨率；
    srcNodata：原来数据的Nodata值
    dstNodata：输出数据的Nodata值
    dstSRS：输出的投影坐标系，可以读取影像的：
            Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
            Projection = Raster.GetProjectionRef()
    resampleAlg:重采样方式,算法包括：
            import gdalconst
            gdalconst.GRA_NearestNeighbour：near
            gdalconst.GRA_Bilinear:bilinear
            gdalconst.GRA_Cubic:cubic
            gdalconst.GRA_CubicSpline:cubicspline
            gdalconst.GRA_Lanczos:lanczos
            gdalconst.GRA_Average:average
            gdalconst.GRA_Mode:mode
    )

"""
import numpy as np
import os,datetime,gdal,gdalconst

# Open datasets，此处的tif文件为模板文件，其提供目标投影坐标
InputImage = r'E:\cq\0-mb\mb.tif'
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()
# 需要重投影的tif文件所在文件夹
inDir = r"F:\重庆数据\潼南-中分影像\32645"+str(i)
# 投影后tif文件存放文件夹
outDir = r'F:\重庆数据\潼南-中分影像\cty/32645'+str(i)

# 如果没有输出文件夹，创建文件夹
if not os.path.exists(outDir):
        os.makedirs(outDir)
inList = [name for name in os.listdir(inDir)]
print(inList)

for file in inList:
    if file.endswith('.tif'):
        input_raster = os.path.join(inDir,file)
    # 输出文件的完整路径
        output_raster = os.path.join(outDir,file)
        if not os.path.exists(output_raster):
            OutTile = gdal.Warp(output_raster, input_raster,
                    dstSRS=Projection,
                    resampleAlg=gdalconst.GRA_NearestNeighbour,
                    dstNodata=0,
                    xRes=10,
                    yRes=10
                    )
        OutTile =None
        print(inList.index(file))
print("END!")

三、矢量样本转栅格

# 样本矢量转栅格
import sys,os
import geopandas as gpd
import rasterio
from rasterio import features
from rasterio.crs import CRS
from rasterio.transform import Affine

class Polygon2Raster:
    def polygon_to_raster(self,shp,raster,pixel,field,field1,code=4326):
        '''
        矢量转栅格
        :param shp: 输入矢量全路径，字符串，无默认值
        :param raster: 输出栅格全路径，字符串，无默认值
        :param pixel: 像元大小，与矢量坐标系相关
        :param field: 栅格像元值字段
        :param Code: 输出坐标系代码，默认为4326
        :return: None
        '''

        # 判断字段是否存在
        shapefile = gpd.read_file(shp)
        if not field in shapefile.columns:
            raise Exception ('输出字段不存在')
        # 判断数据类型
        f_type = shapefile.dtypes.get(field)
        if 'int' in str(f_type):
            shapefile[field] = shapefile[field].astype('int16')
            dtype = 'int16'
        elif 'float' in str(f_type):
            shapefile[field] = shapefile[field].astype('float32')
            dtype = 'float32'
        else:
            raise Exception ('输入字段数据类型为{}，无法进行栅格化操作'.format(f_type))

        bound = shapefile.bounds
        width = int((bound.get('maxx').max()-bound.get('minx').min())/pixel)
        height = int((bound.get('maxy').max()-bound.get('miny').min())/pixel)
        transform = Affine(pixel, 0.0, bound.get('minx').min(),
               0.0, -pixel, bound.get('maxy').max())
        InputImage = r'E:\cq\0-mb\mb.tif'
        Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
        Projection = Raster.GetProjectionRef()
        

        meta = {'driver': 'GTiff',
                'dtype': dtype,
                'nodata': 0,
                'width': width,
                'height': height,
                'count': 2,
                'crs': Raster.GetProjectionRef(),
                'transform': transform}
        print(meta)

        with rasterio.open(raster, 'w+', **meta) as out:
            out_arr = out.read(1)
            shapes = ((geom,value) for geom, value in zip(shapefile.get('geometry'), shapefile.get(field)))
            burned = features.rasterize(shapes=shapes, fill=0, out=out_arr, transform=out.transform)
            out.write_band(1, burned)
            out_arr2 = out.read(2)
            shapes1 = ((geom,value) for geom, value in zip(shapefile.get('geometry'), shapefile.get(field1)))
            burned2 = features.rasterize(shapes=shapes1, fill=0, out=out_arr2, transform=out.transform)
            out.write_band(2, burned2)

if __name__ == '__main__':
    filepath=r"E:\cq\1-yyb"
    ls=os.listdir(filepath)
    for file in ls:
        if file.endswith('.shp'):
            in_shp=str(os.path.join(filepath,file))
            out_tif = r'E:\cq\1-ybbsg/'+os.path.basename(in_shp).strip('.shp')+'.tif'
            pixel_size = 10
            field = 'typecode'
            field1= 'Id'
            Polygon2Raster().polygon_to_raster(in_shp,out_tif,pixel_size,field,field1)
            print(file+'finish')

4. 样本网格扩大到与原图相同的大小

# 样本栅格扩充
i=490313
import gdal
import os
filepath=r'F:\重庆数据\潼南-中分影像\cty\32645'+str(i)
rasterfile=str(os.path.join(filepath,os.listdir(filepath)[0]))
print(rasterfile)
ybfile=r'E:\cq\1-ybbsg\ly2.tif'
out_file=r'E:\cq\1-ybkc/'+str(i)+'.tif'
ds=gdal.Open(rasterfile)
ds2=gdal.Open(ybfile)
if ds2 is None:
    print("打开失败！")
else:
    int_data=ds.GetRasterBand(1)
    yb_data = ds2.GetRasterBand(1).ReadAsArray()
    yb_data1 = ds2.GetRasterBand(2).ReadAsArray()
#     yb_data[yb_data <0]=0
#     yb_data[yb_data==np.nan]=0
    gtif_driver=gdal.GetDriverByName("GTiff")
    out_ds=gtif_driver.Create(out_file,int_data.XSize,int_data.YSize,2,int_data.DataType)
    in_datatf=ds.GetGeoTransform()
    print("moban",in_datatf)
    out_ds.SetProjection(ds.GetProjection())
    out_ds.SetGeoTransform(in_datatf)    
    out_band=out_ds.GetRasterBand(1)
    out_band1=out_ds.GetRasterBand(2)
#     zeros=np.zeros(int_data.XSize,int_data.YSize)
#     out_band.WriteArray(zeros)
    yb_tf=ds2.GetGeoTransform()
    print("yb",yb_tf)
    yb_lx=yb_tf[0]
    yb_ly=yb_tf[3]
    l_x=int(-(in_datatf[0]-yb_lx)/in_datatf[1])
    l_y=int((in_datatf[3]-yb_ly)/-yb_tf[5])
    print(l_x,l_y)
    out_band.WriteArray(yb_data,l_x,l_y)
    out_band1.WriteArray(yb_data1,l_x,l_y)
    print(out_band.XSize,out_band.YSize)
#     out_band[out_band <0]=0
#     out_band[out_data==np.nan]=0

    out_ds.FlushCache()
    del out_ds
    del ds
    del ds2
    print("End!")