0

When using the OGR library or GDAL library with Python script, is it possible to increase the extent of a vector layer without actually adding new data points? In my specific case, I would like to increase the extent of vector layers associated with gpx files so that when I convert them to rasters they all have the same pixel matrix.

EDIT: An attempt of mine to use gdal.Rasterize does not produce a "tiff" file, nor does it cause an error to be reported:

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
4

2 回答 2

1

You need to pass 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) instead of passing the individual kwargs directly. Otherwise, they will be ignored, and the command won't do what you intend. See http://www.gdal.org/python/osgeo.gdal-module.html#RasterizeOptions and http://www.gdal.org/python/ for details and links to the source code (often useful given the terse documentation).

于 2017-02-16T03:31:16.997 回答
0

I was unable to find a method to change the extent of the vector layer. However, I was able to write a python Function that uses gdal.RasterizeLayer() to produce a raster with an extent much larger than the original vector layer. The code for this function is:

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
于 2016-06-28T01:30:43.210 回答