4

I'm trying to create a single raster from several shapefiles. I think what I need to do is to merge many layers into a single 'MEM' object. Then re-project this object to a different coordinate system. The create a raster GeoTIFF with set coordinates and resolutions that will be populated from the 'MEM' object built from my shapefiles.

The problem is that the merging of OGR layers isn't documented very well for the Python module. Do I get the features or layers themselves?

Also this re-projection documentation assumes I want to physically write new shapefiles to disk. I don't see how I would use that for a single merged vector layer.

In my script I have come so far that I populate a list of layers. I have no clue how to merge them together and how to rasterize them. I may have to create an ogr geometry object and rasterize that but I don't know what the "proper" workflow for this is.

UPDATE:

I have made some progress. I have successfully converted shapefiles to new coordinates and outputted them. But when I attempt to rasterize I get this error:

TypeError: in method 'RasterizeLayer', argument 1 of type 'GDALDatasetShadow *'

EDIT: I have updated the save location to be a dataset but the same error persists.

Updated script below:

#Get all files in folder
inputFolderPath = "C:\\input\\"
outputFolderPath = "C:\\output\\"

src_files = glob.glob(inputFolderPath+"\\*.shp")

#Traverse files to export and rasterize
for file in src_files:
    #Convert srs
    unconvertedSHP = gdal.OpenEx(file) # path to input .shp file
    convertedSHP = outputFolderPath+os.path.basename(file)  # path to output .shp file

    ds = gdal.VectorTranslate(convertedSHP, unconvertedSHP, format='ESRI Shapefile', dstSRS='EPSG:3879')

    unconvertedSHP = None
    ds = None

    #Open converted shp files to rasterize
    shape_datasource = ogr.Open(convertedSHP)
    shape_layer = shape_datasource.GetLayer()

    writeTIF = outputFolderPath+ os.path.splitext(os.path.basename(file))[0]+'.tif' #filename with new tif extension

    tiff_ds = gdal.GetDriverByName('GTiff')
    tiff_ds.Create(writeTIF,1024,1024,1,gdal.GDT_UInt16,options = ['PROFILE=GeoTIFF'])     

    #Rasterize    
    ds = gdal.RasterizeLayer(tiff_ds,[1],shape_layer,burn_values=[200], options = ["ALL_TOUCHED=TRUE"]) #The error occurs here

    shape_datasource = None
    ds = None
    unconvertedSHP = None
    convertedSHP = None 
  • According to help(gdal.RasterizeLayer) the first argument (the variable writeTIF in your code) should be a dataset rather than a string. So I recommend you try using gdal.driver.Create to create a blank output first. – lpd2 Aug 4 '17 at 12:17
  • Thank you. I have updated my question but the same error still remains. – Školstvo Aug 4 '17 at 13:50
1

If you have QGIS installed or OSGeo4W.bat, here is one option to merge shps in a directory and re-project them (e.g. reproject from EPSG:27700 to EPSG:4326):

import subprocess

command = 'ogr2ogr -append -update -s_srs EPSG:27700 -t_srs EPSG:4326 -f "ESRI Shapefile" {} {} -nln merge'.format('C:/Temp/mergeDir', 'C:/Temp/shpDir')
input_command = [command]
command_run = subprocess.call(['C:/Program Files/QGIS 2.18/OSGeo4W.bat', input_command], shell=True)

For converting vector to raster see:

http://www.gdal.org/gdal_rasterize.html

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.