Python：用于剪切带有多边形的光栅图像的函数的代码检查

概述我正在请求审核以下代码.我有一个空间参考图像和一个多边形.我编写了一个代码(见下文)来剪辑此图像以保存新图像(剪裁区域).此功能基于要素类的几何图形剪切栅格.基于几何体的剪切意味着您将使用要素类中所有要素的边界来剪切栅格,而不是这些要素的最小边界矩形 输入：多边形图层和一个或多个栅格图层 输出：新的栅格图层,剪切为多边形边界 import osgeo.gdalimport shapefilei 我正在请求审核以下代码.我有一个空间参考图像和一个多边形.我编写了一个代码(见下文)来剪辑此图像以保存新图像(剪裁区域).此功能基于要素类的几何图形剪切栅格.基于几何体的剪切意味着您将使用要素类中所有要素的边界来剪切栅格,而不是这些要素的最小边界矩形

输入：多边形图层和一个或多个栅格图层
输出：新的栅格图层,剪切为多边形边界

import osgeo.gdalimport shapefileimport struct,numpy,pylabimport numpy as npimport ogrimport osr,gdalfrom shapely.geometry import polygonimport osgeo.gdal as gdalimport sysfrom osgeo import gdal,gdalnumeric,ogr,osrimport Image,ImageDrawdef 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    # pan@R_235_4048@    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 shapefilemyshp <- shapefile("shapefile.shp")## Reading the raster you want to cropmyraster <- raster('image.tif')## create a layer only for the shape,the parameter inverse = TRUE or FALSE is impnew_raster = mask(myraster,myshp,filename = "newras.tif",inverse = FALSE)

希望这可以帮助.

总结

以上是内存溢出为你收集整理的Python：用于剪切带有多边形的光栅图像的函数的代码检查全部内容，希望文章能够帮你解决Python：用于剪切带有多边形的光栅图像的函数的代码检查所遇到的程序开发问题。

如果觉得内存溢出网站内容还不错，欢迎将内存溢出网站推荐给程序员好友。