#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
Author:zhaoss
Email:zhaoshaoshuai@hnnydsj.com
Create date:  
File: .py
Description:


Parameters


"""
import os
import glob
import time
import sys
import math
import numpy as np

from osgeo import gdal, ogr, gdalconst

try:
    progress = gdal.TermProgress_nocb
except:
    progress = gdal.TermProgress


def shp2raster(raster_ds, shp_layer, ext):
    # 将行列整数浮点化
    ext = np.array(ext) * 1.0
    # 获取栅格数据的基本信息
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    # 根据最小重叠矩形的范围进行矢量栅格化
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    # 创建mask
    # out = r"F:\test_data\clipraster\gdal_mask2\test3.tif"
    # mask_ds = gdal.GetDriverByName('GTiff').Create(out, int(x_size), int(y_size), 1, gdal.GDT_Byte)
    mask_ds = gdal.GetDriverByName('MEM').Create('', int(x_size), int(y_size), 1, gdal.GDT_Byte)
    mask_ds.SetProjection(raster_prj)
    mask_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    mask_ds.SetGeoTransform(mask_geo)
    # 矢量栅格化
    gdal.RasterizeLayer(mask_ds, [1], shp_layer, burn_values=[1], callback=progress)

    return mask_ds


def min_rect(raster_ds, shp_layer):
    # 获取栅格的大小
    x_size = raster_ds.RasterXSize
    y_size = raster_ds.RasterYSize
    # 获取是矢量的范围
    extent = shp_layer.GetExtent()
    # 获取栅格的放射变换参数
    raster_geo = raster_ds.GetGeoTransform()
    # 计算逆放射变换系数
    raster_inv_geo = gdal.InvGeoTransform(raster_geo)
    # 计算在raster上的行列号
    # 左上
    off_ulx, off_uly = map(round, gdal.ApplyGeoTransform(raster_inv_geo, extent[0], extent[3]))
    # 右下
    off_drx, off_dry = map(round, gdal.ApplyGeoTransform(raster_inv_geo, extent[1], extent[2]))
    # 判断是否有重叠区域
    if off_ulx >= x_size or off_uly >= y_size or off_drx <= 0 or off_dry <= 0:
        sys.exit("Have no overlap")
    # 限定重叠范围在栅格影像上
    # 列
    offset_column = np.array([off_ulx, off_drx])
    offset_column = np.maximum((np.minimum(offset_column, x_size - 1)), 0)
    # 行
    offset_line = np.array([off_uly, off_dry])
    offset_line = np.maximum((np.minimum(offset_line, y_size - 1)), 0)

    return [offset_column[0], offset_line[0], offset_column[1], offset_line[1]]


def mask_raster(raster_ds, mask_ds, outfile, ext):
    # 将行列整数浮点化
    ext = np.array(ext) * 1.0
    # 获取栅格数据的基本信息
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    # 根据最小重叠矩形的范围进行矢量栅格化
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    # 创建输出影像
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size), int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    # 获取掩模
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    # 对原始影像进行掩模并输出
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)

    return 1


def main(raster, shp, out):
    # 打开栅格和矢量影像
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    # 计算矢量和栅格的最小重叠矩形
    offset = min_rect(raster_ds, shp_l)
    # 矢量栅格化
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    # 进行裁剪
    res = mask_raster(raster_ds, mask_ds, out, offset)

    return None


if __name__ == '__main__':
    # 注册所有gdal的驱动
    gdal.AllRegister()
    gdal.SetConfigOption("gdal_FILENAME_IS_UTF8", "YES")
    start_time = time.clock()
    # in_file = r"\\192.168.0.234\nydsj\user\ZSS\gongyi20180416\planet_20180416_gongyi2.tif"
    # shpfile = r"\\192.168.0.234\nydsj\user\ZSS\巩义结果修改\4.dlgq\gongyi_0120.shp"
    # outfile = r"F:\test_data\clipraster\out\gongyi.tif"
    in_file = sys.argv[1]
    shpfile = sys.argv[2]
    outfile = sys.argv[3]
    print(in_file, shpfile, outfile)

    main(in_file, shpfile, outfile)

    end_time = time.clock()

    print("time: %.4f secs." % (end_time - start_time))

01c393de5b7ae7c233205513355566e02134dbde
<mask token>


def mask_raster(raster_ds, mask_ds, outfile, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size),
        int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]
            ), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)
    return 1


def main(raster, shp, out):
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    offset = min_rect(raster_ds, shp_l)
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    res = mask_raster(raster_ds, mask_ds, out, offset)
    return None


<mask token>

<mask token>


def min_rect(raster_ds, shp_layer):
    x_size = raster_ds.RasterXSize
    y_size = raster_ds.RasterYSize
    extent = shp_layer.GetExtent()
    raster_geo = raster_ds.GetGeoTransform()
    raster_inv_geo = gdal.InvGeoTransform(raster_geo)
    off_ulx, off_uly = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[0], extent[3]))
    off_drx, off_dry = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[1], extent[2]))
    if off_ulx >= x_size or off_uly >= y_size or off_drx <= 0 or off_dry <= 0:
        sys.exit('Have no overlap')
    offset_column = np.array([off_ulx, off_drx])
    offset_column = np.maximum(np.minimum(offset_column, x_size - 1), 0)
    offset_line = np.array([off_uly, off_dry])
    offset_line = np.maximum(np.minimum(offset_line, y_size - 1), 0)
    return [offset_column[0], offset_line[0], offset_column[1], offset_line[1]]


def mask_raster(raster_ds, mask_ds, outfile, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size),
        int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]
            ), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)
    return 1


def main(raster, shp, out):
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    offset = min_rect(raster_ds, shp_l)
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    res = mask_raster(raster_ds, mask_ds, out, offset)
    return None


<mask token>

<mask token>


def shp2raster(raster_ds, shp_layer, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    mask_ds = gdal.GetDriverByName('MEM').Create('', int(x_size), int(
        y_size), 1, gdal.GDT_Byte)
    mask_ds.SetProjection(raster_prj)
    mask_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    mask_ds.SetGeoTransform(mask_geo)
    gdal.RasterizeLayer(mask_ds, [1], shp_layer, burn_values=[1], callback=
        progress)
    return mask_ds


def min_rect(raster_ds, shp_layer):
    x_size = raster_ds.RasterXSize
    y_size = raster_ds.RasterYSize
    extent = shp_layer.GetExtent()
    raster_geo = raster_ds.GetGeoTransform()
    raster_inv_geo = gdal.InvGeoTransform(raster_geo)
    off_ulx, off_uly = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[0], extent[3]))
    off_drx, off_dry = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[1], extent[2]))
    if off_ulx >= x_size or off_uly >= y_size or off_drx <= 0 or off_dry <= 0:
        sys.exit('Have no overlap')
    offset_column = np.array([off_ulx, off_drx])
    offset_column = np.maximum(np.minimum(offset_column, x_size - 1), 0)
    offset_line = np.array([off_uly, off_dry])
    offset_line = np.maximum(np.minimum(offset_line, y_size - 1), 0)
    return [offset_column[0], offset_line[0], offset_column[1], offset_line[1]]


def mask_raster(raster_ds, mask_ds, outfile, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size),
        int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]
            ), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)
    return 1


def main(raster, shp, out):
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    offset = min_rect(raster_ds, shp_l)
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    res = mask_raster(raster_ds, mask_ds, out, offset)
    return None


<mask token>

<mask token>
import os
import glob
import time
import sys
import math
import numpy as np
from osgeo import gdal, ogr, gdalconst
try:
    progress = gdal.TermProgress_nocb
except:
    progress = gdal.TermProgress


def shp2raster(raster_ds, shp_layer, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    mask_ds = gdal.GetDriverByName('MEM').Create('', int(x_size), int(
        y_size), 1, gdal.GDT_Byte)
    mask_ds.SetProjection(raster_prj)
    mask_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    mask_ds.SetGeoTransform(mask_geo)
    gdal.RasterizeLayer(mask_ds, [1], shp_layer, burn_values=[1], callback=
        progress)
    return mask_ds


def min_rect(raster_ds, shp_layer):
    x_size = raster_ds.RasterXSize
    y_size = raster_ds.RasterYSize
    extent = shp_layer.GetExtent()
    raster_geo = raster_ds.GetGeoTransform()
    raster_inv_geo = gdal.InvGeoTransform(raster_geo)
    off_ulx, off_uly = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[0], extent[3]))
    off_drx, off_dry = map(round, gdal.ApplyGeoTransform(raster_inv_geo,
        extent[1], extent[2]))
    if off_ulx >= x_size or off_uly >= y_size or off_drx <= 0 or off_dry <= 0:
        sys.exit('Have no overlap')
    offset_column = np.array([off_ulx, off_drx])
    offset_column = np.maximum(np.minimum(offset_column, x_size - 1), 0)
    offset_line = np.array([off_uly, off_dry])
    offset_line = np.maximum(np.minimum(offset_line, y_size - 1), 0)
    return [offset_column[0], offset_line[0], offset_column[1], offset_line[1]]


def mask_raster(raster_ds, mask_ds, outfile, ext):
    ext = np.array(ext) * 1.0
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size),
        int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]
            ), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)
    return 1


def main(raster, shp, out):
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    offset = min_rect(raster_ds, shp_l)
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    res = mask_raster(raster_ds, mask_ds, out, offset)
    return None


if __name__ == '__main__':
    gdal.AllRegister()
    gdal.SetConfigOption('gdal_FILENAME_IS_UTF8', 'YES')
    start_time = time.clock()
    in_file = sys.argv[1]
    shpfile = sys.argv[2]
    outfile = sys.argv[3]
    print(in_file, shpfile, outfile)
    main(in_file, shpfile, outfile)
    end_time = time.clock()
    print('time: %.4f secs.' % (end_time - start_time))

#!/usr/bin/env python
# -*- coding:utf-8 -*-
"""
Author:zhaoss
Email:zhaoshaoshuai@hnnydsj.com
Create date:  
File: .py
Description:


Parameters


"""
import os
import glob
import time
import sys
import math
import numpy as np

from osgeo import gdal, ogr, gdalconst

try:
    progress = gdal.TermProgress_nocb
except:
    progress = gdal.TermProgress


def shp2raster(raster_ds, shp_layer, ext):
    # 将行列整数浮点化
    ext = np.array(ext) * 1.0
    # 获取栅格数据的基本信息
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    # 根据最小重叠矩形的范围进行矢量栅格化
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    # 创建mask
    # out = r"F:\test_data\clipraster\gdal_mask2\test3.tif"
    # mask_ds = gdal.GetDriverByName('GTiff').Create(out, int(x_size), int(y_size), 1, gdal.GDT_Byte)
    mask_ds = gdal.GetDriverByName('MEM').Create('', int(x_size), int(y_size), 1, gdal.GDT_Byte)
    mask_ds.SetProjection(raster_prj)
    mask_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    mask_ds.SetGeoTransform(mask_geo)
    # 矢量栅格化
    gdal.RasterizeLayer(mask_ds, [1], shp_layer, burn_values=[1], callback=progress)

    return mask_ds


def min_rect(raster_ds, shp_layer):
    # 获取栅格的大小
    x_size = raster_ds.RasterXSize
    y_size = raster_ds.RasterYSize
    # 获取是矢量的范围
    extent = shp_layer.GetExtent()
    # 获取栅格的放射变换参数
    raster_geo = raster_ds.GetGeoTransform()
    # 计算逆放射变换系数
    raster_inv_geo = gdal.InvGeoTransform(raster_geo)
    # 计算在raster上的行列号
    # 左上
    off_ulx, off_uly = map(round, gdal.ApplyGeoTransform(raster_inv_geo, extent[0], extent[3]))
    # 右下
    off_drx, off_dry = map(round, gdal.ApplyGeoTransform(raster_inv_geo, extent[1], extent[2]))
    # 判断是否有重叠区域
    if off_ulx >= x_size or off_uly >= y_size or off_drx <= 0 or off_dry <= 0:
        sys.exit("Have no overlap")
    # 限定重叠范围在栅格影像上
    # 列
    offset_column = np.array([off_ulx, off_drx])
    offset_column = np.maximum((np.minimum(offset_column, x_size - 1)), 0)
    # 行
    offset_line = np.array([off_uly, off_dry])
    offset_line = np.maximum((np.minimum(offset_line, y_size - 1)), 0)

    return [offset_column[0], offset_line[0], offset_column[1], offset_line[1]]


def mask_raster(raster_ds, mask_ds, outfile, ext):
    # 将行列整数浮点化
    ext = np.array(ext) * 1.0
    # 获取栅格数据的基本信息
    raster_prj = raster_ds.GetProjection()
    raster_geo = raster_ds.GetGeoTransform()
    bandCount = raster_ds.RasterCount
    dataType = raster_ds.GetRasterBand(1).DataType
    # 根据最小重叠矩形的范围进行矢量栅格化
    ulx, uly = gdal.ApplyGeoTransform(raster_geo, ext[0], ext[1])
    x_size = ext[2] - ext[0]
    y_size = ext[3] - ext[1]
    # 创建输出影像
    result_ds = gdal.GetDriverByName('GTiff').Create(outfile, int(x_size), int(y_size), bandCount, dataType)
    result_ds.SetProjection(raster_prj)
    result_geo = [ulx, raster_geo[1], 0, uly, 0, raster_geo[5]]
    result_ds.SetGeoTransform(result_geo)
    # 获取掩模
    mask = mask_ds.GetRasterBand(1).ReadAsArray()
    mask = 1 - mask
    # 对原始影像进行掩模并输出
    for band in range(bandCount):
        banddata = raster_ds.GetRasterBand(band + 1).ReadAsArray(int(ext[0]), int(ext[1]), int(x_size), int(y_size))
        banddata = np.choose(mask, (banddata, 0))
        result_ds.GetRasterBand(band + 1).WriteArray(banddata)

    return 1


def main(raster, shp, out):
    # 打开栅格和矢量影像
    raster_ds = gdal.Open(raster)
    shp_ds = ogr.Open(shp)
    shp_l = shp_ds.GetLayer()
    # 计算矢量和栅格的最小重叠矩形
    offset = min_rect(raster_ds, shp_l)
    # 矢量栅格化
    mask_ds = shp2raster(raster_ds, shp_l, offset)
    # 进行裁剪
    res = mask_raster(raster_ds, mask_ds, out, offset)

    return None


if __name__ == '__main__':
    # 注册所有gdal的驱动
    gdal.AllRegister()
    gdal.SetConfigOption("gdal_FILENAME_IS_UTF8", "YES")
    start_time = time.clock()
    # in_file = r"\\192.168.0.234\nydsj\user\ZSS\gongyi20180416\planet_20180416_gongyi2.tif"
    # shpfile = r"\\192.168.0.234\nydsj\user\ZSS\巩义结果修改\4.dlgq\gongyi_0120.shp"
    # outfile = r"F:\test_data\clipraster\out\gongyi.tif"
    in_file = sys.argv[1]
    shpfile = sys.argv[2]
    outfile = sys.argv[3]
    print(in_file, shpfile, outfile)

    main(in_file, shpfile, outfile)

    end_time = time.clock()

    print("time: %.4f secs." % (end_time - start_time))
