Geotiff 批量范围裁剪

Geotiff 批量范围裁剪

需求

需要将shpfile中每个多边形范围矩形框的影像裁剪成单独的影像geotiff 格式，影像大小400G+ ，用arcgis对影像压缩后也有10G左右

处理方式

前提 shapefile和影像采用相同的投影坐标系

• a 通过arcpy实现
  1.通过arcgis对多边形进行split。
  2 批量读取shapefile头文件范围,调用arcpy.Clip_management 进行裁剪
  缺点 arcgis分割的shapeflie容易头文件错误导致没法读取正确的范围，从而不能正确完成剩下裁剪工作

shapefile分割

import struct,os
import arcpy

shapeAll="F:\\test\\test\\test.shp"
imageAll="F:\\test\\test"
shapeOutPath="F:\\test\\test\\JG1"
rasterOutPath="F:\\test\\test\\JG2"

if not os.path.exists(shapeOutPath):
    os.makedirs(shapeOutPath)
    print 'create %s ' % shapeOutPath
# 分割shapefile
print 'split shapefile ...'
arcpy.Split_analysis(in_features=shapeAll, split_features=shapeAll, split_field="ID", out_workspace=shapeOutPath, cluster_tolerance="")
print 'split shapefile finished'

批量裁剪

import struct,os
import arcpy

shapeAll="F:\\test\\test\\test.shp"
imageAll="F:\\test\\test"
shapeOutPath="F:\\test\\test\\JG1"
rasterOutPath="F:\\test\\test\\JG2"

if not os.path.exists(shapeOutPath):
    os.makedirs(shapeOutPath)
    print 'create %s ' % shapeOutPath

# get all shapefile
fList = os.listdir(shapeOutPath)
total=0
sucess=0
failed=0
for child in fList:
    if child.endswith(".shp"):
        chunName=os.path.splitext(child)[0]
        imageFile=imageAll +'\\'+'test'+'.tif'
        shapeFile=shapeOutPath + '\\'+child

        # extent
        # print shapeFile
        shp = open(shapeFile)
        header = (struct.Struct(">IIIIIII"),struct.Struct("<IIdddddddd"))
        # SHP   header   ShapFile
        # 0-3   int32  FileCode (0x0000270a)
        # 4-23  int32  Unused 5-uinit32
        # 24-27 int32  FileLength 16bit
        filecode,_,_,_,_,_,shp_file_length = header[0].unpack_from(shp.read(28))
        # 28-31 init32 Version
        # 32-35 init32 ShapeType
        # 36-67 double minX,minY,maxX,maxY
        # 68-83 double minZ,maxZ
        # 84-99 double minM,maxM
        version,_type,lox,loy,hix,hiy,loz,hiz,lom,him = header[1].unpack_from(shp.read(72))

        extent ="%f %f %f %f" % (lox,loy,hix,hiy)
        total=total+1
        if lox<180:
            failed=failed+1
            print ('failed %s ' % shapeFile)
            continue

        # print "start clip iamge : %s  feature : % s" % (imageFile,shapeFile)
        # print "arcpy.Clip_management(in_raster=%s, rectangle=%s, out_raster=%s, in_template_dataset=%s, nodata_value='255', clipping_geometry='NONE', maintain_clipping_extent='NO_MAINTAIN_EXTENT')" % (imageFile,extent,rasterOutPath +'\\'+chunName,shapeFile)
        
        arcpy.Clip_management(in_raster=imageFile, rectangle=extent, out_raster=rasterOutPath +'\\'+chunName+'.tif', in_template_dataset=shapeFile, nodata_value="255", clipping_geometry="NONE", maintain_clipping_extent="NO_MAINTAIN_EXTENT")
        # print "start clip iamge : %s" % shapeFile
        sucess=sucess+1
        #break
print "toatl %d sucess %d failed %d" % (total,sucess,failed)

• b.通过arcpy+postgis 实现
  1 通过qgis或者postgis导入工具把shapefile导入到postgis数据库

-- 根据图斑范围批量裁剪影像 arcpy
drop table if exists sql_tb;
create table sql_tb as 
	select
'arcpy.Clip_management(in_raster="F:/test/test/test.tif", rectangle="'
||st_xmin(geom)||' '||st_ymin(geom)||' '|| st_xmax(geom)||' '||st_ymax(geom)||
'", out_raster="F:/test/test/JG2/'|| id ||'.tif", in_template_dataset="F:/test/test/JG1/'|| id ||'.shp", nodata_value="255",clipping_geometry="NONE", maintain_clipping_extent="NO_MAINTAIN_EXTENT")' as g,'1' as tx
from test
select array_to_string(array_agg(t.g), '')  from sql_tb t group by t.tx 

2 编辑sql 生成Clip_management 的批量处理脚本脚本

arcpy.Clip_management(in_raster="F:/test/test/test.tif", rectangle="107.934060561506 29.4070985335127 107.934285001123 29.4074193314034", out_raster="F:/test/test/JG2/7a03324dc9a644f3abda8c37bdf09ff5.tif", in_template_dataset="F:/test/test/JG1/7a03324dc9a644f3abda8c37bdf09ff5.shp", nodata_value="255", clipping_geometry="NONE")
arcpy.Clip_management(in_raster="F:/test/test/test.tif", rectangle="107.784535772868 29.2162284939819 107.784752923213 29.2165892283663", out_raster="F:/test/test/JG2/c293c44e5662444e8e2225bdc6b0f850.tif", in_template_dataset="F:/test/test/JG1/c293c44e5662444e8e2225bdc6b0f850.shp", nodata_value="255", clipping_geometry="NONE")

c.通过GDAL库读取shapefile中每个多边形范围(或者采用a方案中已经分割的shape)，通过GDAL直接对影像进行裁剪
由于要编译GDAL，切第二种方案基本解决了问题 这种方式没有经过测试

 from osgeo import gdal, gdalnumeric, ogr,struct,os
from PIL import Image, ImageDraw

# 将图像转为数组
def imageToArray(i):
    a = gdalnumeric.fromstring(i.tobytes(), 'b')
    a.shape = i.im.size[1], i.im.size[0]
    return a

# 将数组转为图像
def arrayToImage(a):
    i = Image.frombytes('L', (a.shape[1], a.shape[0]), a.astype('b'))
    return i

# 将经纬度坐标转为像素坐标
def world2Pixel(geoMatrix, x, y):
    ulX = geoMatrix[0]
    ulY = geoMatrix[3]
    xDist = geoMatrix[1]
    yDist = geoMatrix[5]

    pixel = int((x - ulX) / xDist)
    line = int((ulY - y) / yDist)

    return (abs(pixel), abs(line))


# 第一步：准备好要切的tif图、shp文件及输出的文件名
tif_file = "chengdu-xinglong-lake.tif"  # 原始的tif图
shp_file = "xing/xing.shp"  # 要切的shape file
clipped_file = "clipped.tif" # 最后切出来的tif图名字

def clip(tif_file,shp_file,clipped_file):
	# 第二步：将tif文件读取为numpy.ndarray
	srcArray = gdalnumeric.LoadFile(tif_file)

	# 第三步：读取tif文件中的仿信息
	in_ds = gdal.Open(tif_file)
	geoTrans = in_ds.GetGeoTransform()
	print (geoTrans)

	# 第四步：读取shp文件中的图层从而获取要切取部分的extent
	shapef = ogr.Open(shp_file)
	lyr = shapef.GetLayer()
	minX, maxX, minY, maxY = lyr.GetExtent()
	print (111, minX, maxX, minY, maxY)

	# 第五步：将要切取部分的大地坐标转为像素坐标
	ulX, ulY = world2Pixel(geoTrans, minX, maxY)
	lrX, lrY = world2Pixel(geoTrans, maxX, minY)

	# 计算clip出来图的宽和高
	pxWidth = abs(int(lrX - ulX))
	pxHeight = abs(int(lrY - ulY))

	# 从原图中切出来，但这个时候还是矩形的
	clip = srcArray[:, ulY:lrY, ulX:lrX]

	geoTrans = list(geoTrans)

	# 更新左上角的坐标，相当于更新设置切后的tif图的仿射信息
	geoTrans[0] = minX
	geoTrans[3] = maxY

	# 下面主要是构造mask
	points = []
	pixels = []
	poly = lyr.GetNextFeature()
	geom = poly.GetGeometryRef()
	pts = geom.GetGeometryRef(0)
	for p in range(pts.GetPointCount()):
		points.append((pts.GetX(p), pts.GetY(p)))
		print (str(pts.GetX(p))+","+str(pts.GetY(p)))

	for p in points:
		pixels.append(world2Pixel(geoTrans, p[0], p[1]))

	rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)

	rasterize = ImageDraw.Draw(rasterPoly)
	rasterize.polygon(pixels, 0)

	mask = imageToArray(rasterPoly)

	# 按照多边形来切
	clip = gdalnumeric.choose(mask, (clip, 0))

	# 创建切后图的tif文件
	gtif_driver = gdal.GetDriverByName('GTiff')
	out_ds = gtif_driver.Create(clipped_file, pxWidth, pxHeight, 3, in_ds.GetRasterBand(1).DataType)

	# 设置切出来的tif图的仿射信息和参考系统
	out_ds.SetGeoTransform(geoTrans)
	out_ds.SetProjection(in_ds.GetProjection())

	# 写入目标文件
	out_ds.GetRasterBand(1).WriteArray(clip[0])
	out_ds.GetRasterBand(2).WriteArray(clip[1])
	out_ds.GetRasterBand(3).WriteArray(clip[2])

	# 将缓存写入磁盘
	out_ds.FlushCache()
	print("FlushCache succeed")
	
	
shapeAll="F:\\test\\test\\test.shp"
imageAll="F:\\test\\test"
shapeOutPath="F:\\test\\test\\JG1"
rasterOutPath="F:\\test\\test\\JG2"
# get all shapefile
#fList = os.listdir(shapeOutPath.decode('utf-8').encode('gbk'))
fList = os.listdir(shapeOutPath.decode('utf-8').encode('gbk'))
total=0
sucess=0
failed=0
for child in fList:
    if child.endswith(".shp"):
        chunName=os.path.splitext(child)[0]
        imageFile=imageAll +'\\'+'test'+'.tif'
        shapeFile=shapeOutPath + '\\'+child
		clip(imageFile,shapeFile,rasterOutPath+'\\'+chunName+'.tif')