Python:用于剪切带有多边形的光栅图像的函数的代码检查
发布时间:2020-12-16 21:44:12 所属栏目:Python 来源:网络整理
导读:我正在请求审核以下代码.我有一个空间参考图像和一个多边形.我编写了一个代码(见下文)来剪辑此图像以保存新图像(剪裁区域).此功能基于要素类的几何图形剪切栅格.基于几何体的剪切意味着您将使用要素类中所有要素的边界来剪切栅格,而不是这些要素的最小边界矩
我正在请求审核以下代码.我有一个空间参考图像和一个多边形.我编写了一个代码(见下文)来剪辑此图像以保存新图像(剪裁区域).此功能基于要素类的几何图形剪切栅格.基于几何体的剪切意味着您将使用要素类中所有要素的边界来剪切栅格,而不是这些要素的最小边界矩形
输入:多边形图层和一个或多个栅格图层 import osgeo.gdal import shapefile import struct,numpy,pylab import numpy as np import ogr import osr,gdal from shapely.geometry import Polygon import osgeo.gdal as gdal import sys from osgeo import gdal,gdalnumeric,ogr,osr import Image,ImageDraw def world2Pixel(geoMatrix,x,y): """ Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate the pixel location of a geospatial coordinate (http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html) geoMatrix [0] = top left x (x Origin) [1] = w-e pixel resolution (pixel Width) [2] = rotation,0 if image is "north up" [3] = top left y (y Origin) [4] = rotation,0 if image is "north up" [5] = n-s pixel resolution (pixel Height) """ ulX = geoMatrix[0] ulY = geoMatrix[3] xDist = geoMatrix[1] yDist = geoMatrix[5] rtnX = geoMatrix[2] rtnY = geoMatrix[4] pixel = np.round((x - ulX) / xDist).astype(np.int) line = np.round((ulY - y) / xDist).astype(np.int) return (pixel,line) def Pixel2world(geoMatrix,y): ulX = geoMatrix[0] ulY = geoMatrix[3] xDist = geoMatrix[1] yDist = geoMatrix[5] coorX = (ulX + (x * xDist)) coorY = (ulY + (y * yDist)) return (coorX,coorY) def RASTERClipByPolygon(inFile,poly,outFile): # Open the image as a read only image ds = osgeo.gdal.Open(inFile,gdal.GA_ReadOnly) # Check the ds (=dataset) has been successfully open # otherwise exit the script with an error message. if ds is None: raise SystemExit("The raster could not openned") # Get image georeferencing information. geoMatrix = ds.GetGeoTransform() ulX = geoMatrix[0] ulY = geoMatrix[3] xDist = geoMatrix[1] yDist = geoMatrix[5] rtnX = geoMatrix[2] rtnY = geoMatrix[4] # get the WKT (= Well-known text) dsWKT = ds.GetProjectionRef() # get driver information DriverName = ds.GetDriver().ShortName # open shapefile (= border of are of interest) shp = osgeo.ogr.Open(poly) if len(shp.GetLayer()) != 1: raise SystemExit('The shapefile must have exactly one layer') # Create an OGR layer from a boundary shapefile layer = shp.GetLayer(0) feature = layer.GetNextFeature() geometry = feature.GetGeometryRef() # Make sure that it is a polygon if geometry.GetGeometryType() != osgeo.ogr.wkbPolygon: raise SystemExit('This module can only load polygon') # get Extent of the clip area X_min,X_max,Y_min,Y_max = layer.GetExtent() # Convert the layer extent to image pixel coordinates uldX,uldY = world2Pixel(geoMatrix,X_min,Y_max) lrdX,lrdY = world2Pixel(geoMatrix,Y_min) # Calculate the pixel size of the new image pxWidth = int(lrdX - uldX) pxHeight = int(lrdY - uldY) # get the Coodinate of left-up vertex of the pixel X_minPixel,Y_maxPixel = Pixel2world(geoMatrix,uldX,uldY) # get polygon's vertices pts = geometry.GetGeometryRef(0) points = [] for p in range(pts.GetPointCount()): points.append((pts.GetX(p),pts.GetY(p))) pnts = np.array(points).transpose() # work band by band nBands = ds.RasterCount # panchromatic if nBands == 1: band = ds.GetRasterBand(1) # get nodata value nodata = band.GetNoDataValue() # convert band in Array bandArray = band.ReadAsArray() del band # clip arrey bandArray_Area = bandArray[uldY:lrdY,uldX:lrdX] del bandArray # Create 2D Polygon Mask. Mode 'L',not '1',because # Numpy-1.5.0 / PIL-1.1.7 does not support the numpy.array(img) # conversion nicely for bivalue images. img = Image.new('L',(pxWidth,pxHeight),0) target_ds = gdal.GetDriverByName(DriverName).Create(outFile,pxWidth,pxHeight,nBands,ds.GetRasterBand(1).DataType) target_ds.SetGeoTransform((X_minPixel,xDist,rtnX,Y_maxPixel,rtnY,yDist)) pixels,line = world2Pixel(target_ds.GetGeoTransform(),pnts[0],pnts[1]) listdata = [(pixels[i],line[i]) for i in xrange(len(pixels))] ImageDraw.Draw(img).polygon(listdata,outline=1,fill=1) mask = numpy.array(img) bandArray_Masked = bandArray_Area*mask del bandArray_Area,mask target_ds.GetRasterBand(nBands).WriteArray(bandArray_Masked) target_ds.GetRasterBand(nBands).SetNoDataValue(nodata) else: img = Image.new('L',ds.GetRasterBand(1).DataType) target_ds.SetGeoTransform((X_min,Y_max,fill=1) mask = numpy.array(img) for bandno in range(1,nBands+1): band = ds.GetRasterBand(bandno) nodata = band.GetNoDataValue() # convert band in Array bandArray = band.ReadAsArray() del band # clip arrey bandArray_Area = bandArray[ulY:lrY,ulX:lrX] del bandArray bandArray_Masked = bandArray_Area*mask target_ds.GetRasterBand(bandno).WriteArray(bandArray_Masked) del bandArray_Area target_ds.GetRasterBand(bandno).SetNoDataValue(nodata) # set the reference info if len(dsWKT) is 0: # Source has no projection (needs GDAL >= 1.7.0 to work) target_ds.SetProjection('LOCAL_CS["arbitrary"]') else: # Make the target raster have the same projection as the source target_ds.SetProjection(dsWKT) target_ds = None 解决方法
这可以在R中轻松完成.我刚刚意识到问题是针对python的.因此我做了编辑.有包装可用于在R中的python或python中运行R,请检查包rpy2.
## Read the shapefile myshp <- shapefile("shapefile.shp") ## Reading the raster you want to crop myraster <- raster('image.tif') ## create a layer only for the shape,the parameter inverse = TRUE or FALSE is imp new_raster = mask(myraster,myshp,filename = "newras.tif",inverse = FALSE) 希望这可以帮助. (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |