2016-06-24 16 views
0

Bei Verwendung der OGR-Bibliothek oder GDAL-Bibliothek mit Python-Skript ist es möglich, den Umfang einer Vektorschicht zu erhöhen, ohne dass tatsächlich neue Datenpunkte hinzugefügt werden? In meinem speziellen Fall möchte ich das Ausmaß der Vektorebenen erhöhen, die mit GPX-Dateien verbunden sind, so dass, wenn ich sie in Raster umwandle, sie alle dieselbe Pixelmatrix haben.Vergrößern der Vektorschicht mit OGR oder GDAL?

EDIT: Ein Versuch, von mir gdal.Rasterize verwendet keine „tiff“ Datei erzeugt, noch einen Fehler verursacht werden, berichtet:

import os 
import gdal 
import ogr  
import math 

os.chdir(r'C:\Users\pipi\Documents\Rogaine\Tarlo\gpx') #folder containing gpx files 
vector_fn = '6_hour_Autumngaine_w_Tom_Elle.gpx' #filename of input gpxfile 
pixel_size = 20 #units are in m if gpx file is left in wgs84 
raster_fn = '0011a.tif' # Filename of the raster Tiff that will be created 

driver = ogr.GetDriverByName('GPX') 
source_ds = driver.Open(vector_fn, 0) 
source_layer = source_ds.GetLayer('track_points') #returns the 'track points' layer of the data source 
SR = source_layer.GetSpatialRef().ExportToWkt() 

#_______USING VALUES FROM THE FILE___________ 
x_min1, x_max1, y_min1, y_max1 = source_layer.GetExtent() 

pixel_sizey = pixel_size/(111.2*math.pow(10,3)) #determines an approximate x and y size because of geographic coordinates. 
pixel_sizex = pixel_size/(math.cos(((y_max1 + y_min1)/2)*(math.pi/180))*111.2*math.pow(10,3)) 
print (pixel_sizey, pixel_sizex) 
x_res = int((x_max1 - x_min1)/pixel_sizex) 
y_res = int((y_max1 - y_min1)/pixel_sizey) 
print (x_res, y_res) 

layer_list = ['track_points'] 

gdal.Rasterize(raster_fn, vector_fn, format='GTiff', outputBounds=[x_min1, y_min1, x_max1, y_max1], outputSRS=SR, xRes=x_res, yRes=y_res, burnValues=[1], layers=layer_list) 

target_ds = None 
vector_fn = None 
source_layer = None 
source_ds = None 
+0

Verwenden Sie 'gdal.Rasterize()'? – Benjamin

+0

Nein, ich verwende 'gdal.RasterizeLayer()'. –

+0

Ich würde gerne Hilfe oder eine Kopie von Python-Code, der mit gdal.Rasterize() funktioniert, auch wenn es für einen anderen Zweck ist. Meine Bemühungen, gdal.Rasterize() zu verwenden, führen nicht dazu, dass kein tiff erzeugt wird. –

Antwort

1

Sie müssen options=gdal.RasterizeOptions(format='GTiff', outputBounds=[x_min1, y_min1, x_max1, y_max1], outputSRS=SR, xRes=x_res, yRes=y_res, burnValues=[1], layers=layer_list) passieren, anstatt die einzelnen kwargs direkt zu übergeben. Andernfalls werden sie ignoriert, und der Befehl wird nicht das tun, was Sie beabsichtigen. Einzelheiten und Links zum Quellcode finden Sie unter http://www.gdal.org/python/osgeo.gdal-module.html#RasterizeOptions und http://www.gdal.org/python/ (oft hilfreich in der kurzen Dokumentation).

0

Ich war nicht in der Lage, ein Verfahren zu finden, das Ausmaß zu ändern der Vektorschicht. Ich konnte jedoch eine Python-Funktion schreiben, die gdal.RasterizeLayer() verwendet, um ein Raster mit einer Ausdehnung zu erzeugen, die viel größer ist als die ursprüngliche Vektorebene. Der Code für diese Funktion ist:

import os 
import gdal 
import ogr  

def RasterizeLarge(name, layer, extent, pixel_size): 
    """Used to rasterize a layer where the raster extent is much larger than the layer extent 
    Arguments: 
    name  -- (string) filename without extension of raster to be produced 
    layer  -- (vector layer object) vector layer containing the data to be rasterized (tested with point data) 
    extent  -- (list: x_min, x_max, y_min, y_max) extent of raster to be produced 
    pixel_size -- (list: x_pixel_size, y_pixel_size) 1 or 2 pixel different pixel sizes may be sent 
    """ 

    if isinstance(pixel_size, (list, tuple)): 
     x_pixel_size = pixel_size[0] 
     y_pixel_size = pixel_size[1] 

    else: 
     x_pixel_size = y_pixel_size = pixel_size 

    x_min, x_max, y_min, y_max = extent  
    # determines the x and y resolution of the file (lg = large) 
    x_res_lg = int((x_max - x_min)/x_pixel_size)+2 
    y_res_lg = int((y_max - y_min)/y_pixel_size)+2 

    if x_res_lg > 1 and y_res_lg > 1: 
     pass 
    else: 
     print ('Your pixel size is larger than the extent in one dimension or more') 
     return 

    x_min_sm, x_max_sm, y_min_sm, y_max_sm = layer.GetExtent() 

    if x_min_sm > x_min and x_max_sm < x_max and y_min_sm > y_min and y_max_sm < y_max: 
     pass 
    else: 
     print ('The extent of the layer is in one or more parts outside of the extent provided') 
     return 

    nx = int((x_min_sm - x_min)/x_pixel_size) #(number of pixels between main raster origin and minor raster) 
    ny = int((y_max - y_max_sm)/y_pixel_size) 

    x_res_sm = int((x_max_sm - x_min_sm)/x_pixel_size)+2 
    y_res_sm = int((y_max_sm - y_min_sm)/y_pixel_size)+2 

    #determines upper left corner of small layer raster 
    x_min_sm = x_min + nx * x_pixel_size 
    y_max_sm = y_max - ny * y_pixel_size 

    #______Creates a temporary raster file for the small raster__________ 
    try: 
     # create the target raster file with 1 band 
     sm_ds = gdal.GetDriverByName('GTiff').Create('tempsmall.tif', x_res_sm, y_res_sm, 1, gdal.GDT_Byte) 
     sm_ds.SetGeoTransform((x_min_sm, x_pixel_size, 0, y_max_sm, 0, -y_pixel_size)) 
     sm_ds.SetProjection(layer.GetSpatialRef().ExportToWkt()) 
     gdal.RasterizeLayer(sm_ds, [1], layer, burn_values=[1]) 

     sm_ds.FlushCache() 
     #______Gets data from the new raster in the form of an array________ 
     in_band = sm_ds.GetRasterBand(1) 
     in_band.SetNoDataValue(0) 
     sm_data = in_band.ReadAsArray() 
    finally: 
     sm_ds = None #flushes data from memory. Without this you often get an empty raster. 

    #_____Creates an output file with the provided name and extent that contains the small raster. 
    name = name + '.tif' 

    try: 
     lg_ds = gdal.GetDriverByName('GTiff').Create(name, x_res_lg, y_res_lg, 1, gdal.GDT_Byte) 

     if lg_ds is None: 
      print 'Could not create tif' 
      return 
     else: 
     pass 

     lg_ds.SetProjection(layer.GetSpatialRef().ExportToWkt()) 
     lg_ds.SetGeoTransform((x_min, x_pixel_size, 0.0, y_max, 0.0, -y_pixel_size)) 
     lg_band = lg_ds.GetRasterBand(1) 
     lg_data = in_band.ReadAsArray() 

     lg_band.WriteArray(sm_data, xoff = nx, yoff = ny) 
     lg_band.SetNoDataValue(0) 
     lg_band.FlushCache() 
     lg_band.ComputeStatistics(False) 
     lg_band = None 

    finally: 
     del lg_ds, lg_band, in_band 
     os.remove('tempsmall.tif') 
    return 
Verwandte Themen