python에서 gdal, ogr 등으로 된 shapefile이있는 GTiff 마스크

Question

괜찮습니다. 약간의 재미 만보고 난 후에 두 번째 주석 행의 사이트 하이퍼 링크에서 스크립트를 수정했습니다. 스크립트의 목적은 여러 폴리곤 (각각 "Name"레코드가있는)을 가진 shapefile을 가진 GTiff 형식의 큰 래스터 (32 비트 파이썬 2.7.5 어플리케이션에 적합하지 않음)를 클립/마스크하고, 클리핑 된 래스터를 "클립"하위 디렉토리에 저장합니다. 여기서 각 마스크 된 그리드는 각 다각형의 "이름"다음에 이름이 지정됩니다. 원본 스크립트와 마찬가지로, GTiff와 shapefile은 동일한 프로젝션에 있고 올바르게 겹쳐져 있으며 ~ 100 개의 마스크/초를 처리한다고 가정합니다. 그러나, 나는 그 스크립트를 수정하여 1) float 값의 표고 그리드로 작업하고, 2) 큰 폴리곤의 창을 현재 폴리곤으로 제한된 메모리로로드하고 (즉, 메모리로드를 줄이기 위해), 2) 익스포트 GTiff의 위치 정보와 가치는 맞다.

그러나, 나는 "오른쪽 그림자"라고 부를 것입니다 각 마스크 그리드에 문제가 있습니다. 그것은 그 라인의 오른쪽이 주어진 다각형의 바깥쪽에있는 다각형의 모든 수직선을위한 것이고, 마스킹 된 격자는 그 다각형면의 오른쪽에 하나의 래스터 셀을 포함합니다.

따라서 제 질문은 마스크 된 격자에 오른쪽 그림자를주는 것은 제가 잘못하고있는 것입니까?

다른 모양을 재현 할 수 있도록 예제 shapefile 및 tif를 게시하는 방법을 알아 내려고 노력할 것입니다. 또한 아래 코드에는 정수 값 위성 이미지에 대한 주석 행이 있습니다 (예 : geospatialpython.com의 원본 코드와 동일).

# RasterClipper.py - clip a geospatial image using a shapefile 
# http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html 
# http://gis.stackexchange.com/questions/57005/python-gdal-write-new-raster-using-projection-from-old 

import os, sys, time, Tkinter as Tk, tkFileDialog 
import operator 
from osgeo import gdal, gdalnumeric, ogr, osr 
import Image, ImageDraw 

def SelectFile(req = 'Please select a file:', ft='txt'): 
    """ Customizable file-selection dialogue window, returns list() = [full path, root path, and filename]. """ 
    try: # Try to select a csv dataset 
     foptions = dict(filetypes=[(ft+' file','*.'+ft)], defaultextension='.'+ft) 
     root = Tk.Tk(); root.withdraw(); fname = tkFileDialog.askopenfilename(title=req, **foptions); root.destroy() 
     return [fname]+list(os.path.split(fname)) 
    except: print "Error: {0}".format(sys.exc_info()[1]); time.sleep(5); sys.exit() 

def rnd(v, N): return int(round(v/float(N))*N) 
def rnd2(v): return int(round(v)) 

# Raster image to clip 
rname = SelectFile('Please select your TIF DEM:',ft='tif') 
raster = rname[2] 
print 'DEM:', raster 
os.chdir(rname[1]) 

# Polygon shapefile used to clip 
shp = SelectFile('Please select your shapefile of catchments (requires Name field):',ft='shp')[2] 
print 'shp:', shp 

qs = raw_input('Do you want to stretch the image? (y/n)') 
qs = True if qs == 'y' else False 

# Name of base clip raster file(s) 
if not os.path.exists('./clip/'): os.mkdir('./clip/') 
output = "/clip/clip" 

# This function will convert the rasterized clipper shapefile 
# to a mask for use within GDAL. 
def imageToArray(i): 
    """ 
    Converts a Python Imaging Library array to a 
    gdalnumeric image. 
    """ 
    a=gdalnumeric.fromstring(i.tostring(),'b') 
    a.shape=i.im.size[1], i.im.size[0] 
    return a 

def arrayToImage(a): 
    """ 
    Converts a gdalnumeric array to a 
    Python Imaging Library Image. 
    """ 
    i=Image.fromstring('L',(a.shape[1],a.shape[0]), (a.astype('b')).tostring()) 
    return i 

def world2Pixel(geoMatrix, x, y, N= 5, r=True): 
    """ 
    Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate 
    the pixel location of a geospatial coordinate 
    """ 
    ulX = geoMatrix[0] 
    ulY = geoMatrix[3] 
    xDist = geoMatrix[1] 
    yDist = geoMatrix[5] 
    rtnX = geoMatrix[2] 
    rtnY = geoMatrix[4] 
    if r: 
     pixel = int(round(x - ulX)/xDist) 
     line = int(round(ulY - y)/xDist) 
    else: 
     pixel = int(rnd(x - ulX, N)/xDist) 
     line = int(rnd(ulY - y, N)/xDist) 
    return (pixel, line) 

def histogram(a, bins=range(0,256)): 
    """ 
    Histogram function for multi-dimensional array. 
    a = array 
    bins = 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] 
    return hist 

def stretch(a): 
    """ 
    Performs a histogram stretch on a gdalnumeric array image. 
    """ 
    hist = histogram(a) 
    im = arrayToImage(a) 
    lut = [] 
    for b in range(0, len(hist), 256): 
     # step size 
     step = reduce(operator.add, hist[b:b+256])/255 
     # create equalization lookup table 
     n = 0 
     for i in range(256): 
      lut.append(n/step) 
      n = n + hist[i+b] 
    im = im.point(lut) 
    return imageToArray(im) 

# Also load as a gdal image to get geotransform 
# (world file) info 
srcImage = gdal.Open(raster) 
geoTrans_src = srcImage.GetGeoTransform() 
#print geoTrans_src 
pxs = int(geoTrans_src[1]) 
srcband = srcImage.GetRasterBand(1) 
ndv = -9999.0 
#ndv = 0 

# Create an OGR layer from a boundary shapefile 
shapef = ogr.Open(shp) 
lyr = shapef.GetLayer() 
minXl, maxXl, minYl, maxYl = lyr.GetExtent() 
ulXl, ulYl = world2Pixel(geoTrans_src, minXl, maxYl) 
lrXl, lrYl = world2Pixel(geoTrans_src, maxXl, minYl) 
#poly = lyr.GetNextFeature() 
for poly in lyr: 
    pnm = poly.GetField("Name") 

    # Convert the layer extent to image pixel coordinates 
    geom = poly.GetGeometryRef() 
    #print geom.GetEnvelope() 
    minX, maxX, minY, maxY = geom.GetEnvelope() 

    geoTrans = geoTrans_src 
    ulX, ulY = world2Pixel(geoTrans, minX, maxY) 
    lrX, lrY = world2Pixel(geoTrans, maxX, minY) 

    # Calculate the pixel size of the new image 
    pxWidth = int(lrX - ulX) 
    pxHeight = int(lrY - ulY) 

    # Load the source data as a gdalnumeric array 
    #srcArray = gdalnumeric.LoadFile(raster) 
    clip = gdalnumeric.BandReadAsArray(srcband, xoff=ulX, yoff=ulY, win_xsize=pxWidth, win_ysize=pxHeight) 
    #clip = srcArray[:, ulY:lrY, ulX:lrX] 

    # Create a new geomatrix for the image 
    geoTrans = 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) 
    for p in range(pts.GetPointCount()): 
     points.append((pts.GetX(p), 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 the image using the mask 
    #clip = gdalnumeric.choose(mask, (clip, 0)).astype(gdalnumeric.uint8) 
    clip = gdalnumeric.choose(mask, (clip, ndv)).astype(gdalnumeric.numpy.float) 

    # This image has 3 bands so we stretch each one to make them 
    # visually brighter 
    #for i in range(3): 
    # clip[i,:,:] = stretch(clip[i,:,:]) 
    if qs: clip[:,:] = stretch(clip[:,:]) 

    # Save ndvi as tiff 
    outputi = rname[1]+output+'_'+pnm+'.tif' 
    #gdalnumeric.SaveArray(clip, outputi, format="GTiff", prototype=srcImage) 
    driver = gdal.GetDriverByName('GTiff') 
    DataSet = driver.Create(outputi, pxWidth, pxHeight, 1, gdal.GDT_Float64) 
    #DataSet = driver.Create(outputi, pxWidth, pxHeight, 1, gdal.GDT_Int32) 
    DataSet.SetGeoTransform(geoTrans) 
    Projection = osr.SpatialReference() 
    Projection.ImportFromWkt(srcImage.GetProjectionRef()) 
    DataSet.SetProjection(Projection.ExportToWkt()) 
    # Write the array 
    DataSet.GetRasterBand(1).WriteArray(clip) 
    DataSet.GetRasterBand(1).SetNoDataValue(ndv) 

    # Save ndvi as an 8-bit jpeg for an easy, quick preview 
    #clip = clip.astype(gdalnumeric.uint8) 
    #gdalnumeric.SaveArray(clip, rname[1]+outputi+'.jpg', format="JPEG") 
    #print '\t\tSaved:', outputi, '-.tif, -.jpg' 
    print 'Saved:', outputi 
    del mask, clip, geom 
    del driver, DataSet 

del shapef, srcImage, srcband 

Answer 1

이 기능은 이미 gdal 명령 줄 유틸리티에 통합되어 있습니다. 귀하의 케이스를 감안할 때, 파이썬에서이 작업을 직접하고 싶지는 않습니다.

OGR을 사용하여 기하 구조를 반복하고 각 매개 변수에 대해 gdalwarp을 적절한 매개 변수로 호출 할 수 있습니다.

import ogr 
import subprocess 

inraster = 'NE1_HR_LC_SR_W_DR\NE1_HR_LC_SR_W_DR.tif' 
inshape = '110m_cultural\ne_110m_admin_0_countries_lakes.shp' 

ds = ogr.Open(inshape) 
lyr = ds.GetLayer(0) 

lyr.ResetReading() 
ft = lyr.GetNextFeature() 

while ft: 

    country_name = ft.GetFieldAsString('admin') 

    outraster = inraster.replace('.tif', '_%s.tif' % country_name.replace(' ', '_'))  
    subprocess.call(['gdalwarp', inraster, outraster, '-cutline', inshape, 
        '-crop_to_cutline', '-cwhere', "'admin'='%s'" % country_name]) 

    ft = lyr.GetNextFeature() 

ds = None 

내가 위의 예에서 자연 지구에서 몇 가지 예제 데이터를 사용하고, 브라질의 컷 아웃과 같이 보인다 :

만의 영역으로 이미지를 자르려면 polygon의 마스크가 포함되지 않도록 Shapefile을 변환 할 수 있습니다. 또는 shapefile을 풀고 gdal_translate을 -projwin으로 호출하여 관심 영역을 지정하십시오.