python切割图像_python图像切割

# -*- coding: utf-8 -*-"""@author lzugis@date 2017-06-02@brief 利用shp裁剪影像"""fromosgeo importgdal, gdalnumeric, ogr,gdal_array

fromPIL importImage, ImageDraw

importos

importoperator

fromfunctools importreduce

importnumpy asnp

gdal.UseExceptions()

# This function will convert the rasterized clipper shapefile# to a mask for use within GDAL.defimageToArray(i):

"""Converts a Python Imaging Library array to agdalnumeric image."""a = gdalnumeric.fromstring(i.tobytes(), 'b')

a.shape = i.im.size[1], i.im.size[0]

returna

defarrayToImage(a):

"""Converts a gdalnumeric array to aPython Imaging Library Image."""i = Image.frombytes('L', (a.shape[1], a.shape[0]),

(a.astype('b')).tobytes())

returni

defworld2Pixel(geoMatrix, x, y):

"""Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculatethe pixel location of a geospatial coordinate"""ulX = geoMatrix[0]

ulY = geoMatrix[3]

xDist = geoMatrix[1]

pixel = int((x - ulX) / xDist)

line = int((ulY - y) / xDist)

return(pixel, line)

## EDIT: this is basically an overloaded# version of the gdal_array.OpenArray passing in xoff, yoff explicitly# so we can pass these params off to CopyDatasetInfo#defOpenArray(array, prototype_ds=None, xoff=0, yoff=0):

ds = gdal.Open(gdalnumeric.GetArrayFilename(array))

ifds is not None andprototype_ds is not None:

iftype(prototype_ds).__name__== 'str':

prototype_ds = gdal.Open(prototype_ds)

ifprototype_ds is not None:

gdalnumeric.CopyDatasetInfo(prototype_ds, ds, xoff=xoff, yoff=yoff)

returnds

defhistogram(a, bins=range(0, 256)):

"""Histogram function for multi-dimensional array.a = arraybins = range of numbers to match"""fa = a.flat

n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)

n = gdalnumeric.concatenate([n, [len(fa)]])

hist = n[1:] - n[:-1]

returnhist

defstretch(a):

"""Performs a histogram stretch on a gdalnumeric array image."""real = np.zeros(a.shape)

fori inrange(len(a)):

forj inrange(len(a[i])):

val = a[i][j].real

real[i, j] = val

hist = histogram(real)

im = arrayToImage(real)

lut = []

forb inrange(0, len(hist), 256):

# step sizestep = reduce(operator.add, hist[b:b+256]) / 255# create equalization lookup tablen = 0fori inrange(256):

lut.append(n / step)

n = n + hist[i+b]

im = im.point(lut)

returnimageToArray(im)

importtime

defdraw(points):

img = Image.fromarray(points)

file_path2 = "{path}/ncData/{filename}".format(#注意网站和程序的目录不一样path=sys.path[0],

filename= str(time.time())+'.tif')

img.save(file_path2)

fromscipy importmisc

defmain(shapefile_path, raster_path):

# Load the source data as a gdalnumeric arraysrcArray = gdalnumeric.LoadFile(raster_path)

# Also load as a gdal image to get geotransform# (world file) infosrcImage = gdal.Open(raster_path)

geoTrans = srcImage.GetGeoTransform()

# Create an OGR layer from a boundary shapefileshapef = ogr.Open(shapefile_path)

lyr = shapef.GetLayer(os.path.split(os.path.splitext(shapefile_path)[0])[1])

poly = lyr.GetNextFeature()

# Convert the layer extent to image pixel coordinatesminX, maxX, minY, maxY = lyr.GetExtent()

ulX, ulY = world2Pixel(geoTrans, minX, maxY)

lrX, lrY = world2Pixel(geoTrans, maxX, minY)

# Calculate the pixel size of the new imagepxWidth = int(lrX - ulX)

pxHeight = int(lrY - ulY)

clip = srcArray[ulY:lrY, ulX:lrX]

## EDIT: create pixel offset to pass to new image Projection info#xoffset = ulX

yoffset = ulY

print("Xoffset, Yoffset = ( %f, %f )"% (xoffset, yoffset))

# Create a new geomatrix for the imagegeoTrans = list(geoTrans)

geoTrans[0] = minX

geoTrans[3] = maxY

# Map points to pixels for drawing the# boundary on a blank 8-bit,# black and white, mask image.points = []

pixels = []

geom = poly.GetGeometryRef()

pts = geom.GetGeometryRef(0)

forp inrange(pts.GetPointCount()):

points.append((pts.GetX(p), pts.GetY(p)))

forp inpoints:

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 the image using the mask# clip = gdalnumeric.choose(mask, \# (clip, 0)).astype(gdalnumeric.uint8)clip = gdalnumeric.choose(mask, \

(clip, 0))

real = np.zeros(clip.shape)

fori inrange(len(clip)):

forj inrange(len(clip[i])):

val = clip[i][j].real

if(val < -4):

val = -4real[i, j] = val

dst_ds = gdal.GetDriverByName('GTiff').Create("hello86568.tif", 300, 300, 1, gdal.GDT_CFloat32)

# dst_ds.SetGeoTransform([444720, 30, 0, 3751320, 0, -30])raster= np.zeros(real.shape, dtype=np.float32)

dst_ds.GetRasterBand(1).WriteArray(real)

# Once we're done, close properly the datasetdst_ds= None# This image has 3 bands so we stretch each one to make them# visually brighter# for i in range(3):clip2 = stretch(clip)

# Save new tiff## EDIT: instead of SaveArray, let's break all the# SaveArray steps out more explicity so# we can overwrite the offset of the destination# raster#### the old way using SaveArray## gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)###### gtiffDriver = gdal.GetDriverByName('GTiff')# if gtiffDriver is None:# raise ValueError("Can't find GeoTiff Driver")# gtiffDriver.CreateCopy("beijing9.tif",# gdal_array.OpenArray(clip, prototype_ds=raster_path)# )# real = np.zeros(clip.shape)# for i in range(len(clip)):# for j in range(len(clip[i])):# val = clip[i][j].real# real[i, j] = val# draw(real)# Save as an 8-bit jpeg for an easy, quick previewclip3 = clip2.astype(gdalnumeric.uint8)

gdalnumeric.SaveArray(clip3, "beijing.jpg", format="JPEG")

# misc.imsave("beijing7.png", clip2)gdal.ErrorReset()

importsys

if__name__ == '__main__':

shapefile_path = "{path}/temp/{filename}".format( # 注意网站和程序的目录不一样 且必须是完整的路径path=sys.path[0],

filename='donghai.shp')

raster_path = "{path}/temp/{filename}".format( # 注意网站和程序的目录不一样 且必须是完整的路径path=sys.path[0],

filename='point1522454253.tif')

# shapefile_path, raster_path# shapefile_path = 'temp/donghai.shp'# raster_path = 'temp/point1522361587.tif'main(shapefile_path, raster_path)