Hello everyone, I have another article. Recently, my colleagues were working on ecological service-related projects and needed to operate vector data, so I checked the relevant information. Today I will share with you how to use Python’s GDAL library to implement raster to feature and feature to raster ( Convert raster and vector to each other). In fact, I have shared the code for converting raster to polygon and calculating feature area before. If you are interested, you can check it out: [Python & GIS] GDAL raster to polygon & calculate vector area
1. Raster to surface
Here the value corresponding to the raster is written into the attribute field Value of the vector. The notes are very detailed and don’t say much.
def Raster_To_Vector(input_raster, output_shp):
"""
:param input_raster: Enter the raster data to be converted
:param output_shp: path of output vector
:return:None
"""
print("Converting raster to polygon...")
# Convert raster to polygon
ds_raster = gdal.Open(input_raster) # Read raster data in the path
band_raster = ds_raster.GetRasterBand(1) # Get the band that needs to be converted to vector
proj_raster = osr.SpatialReference()
proj_raster.ImportFromWkt(ds_raster.GetProjection())
# Assign the projection information of raster data to vector
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(output_shp): # If the file already exists, delete it and do it again
driver.DeleteDataSource(output_shp)
polygon = driver.CreateDataSource(output_shp) # Create data resources
layer_polygon = polygon.CreateLayer("Shp", srs=proj_raster, geom_type=ogr.wkbMultiPolygon) # Create a layer and define multiple faces
new_field = ogr.FieldDefn('value', ogr.OFTReal) # Add a field to the target shp file to store the pixel value of the original raster
layer_polygon.CreateField(new_field)
gdal.FPolygonize(band_raster, None, layer_polygon, 0) # Core function, which performs raster to vector operation
polygon.SyncToDisk()
polygon = None
2. Area to Raster
There’s a big problem here! ! ! The parameters in the conversion function in GDAL require reference values, so we need to enter a raster as a reference to convert the vector into a raster. This results in a very limited application of this function. However, when converting an area to a raster in GIS, you do not need to refer to the raster, and you can even choose the raster size.
Speechless…, but this function is not useless. This function can convert our deep learning samples into rasters. Do you understand what I mean? When we draw samples, they are all vectors, but the input of some models requires raster data. At this time, this function is still somewhat useful.
def Vector_To_Raster(input_shp, refer_raster, attribute, output_raster):
"""
:param input_shp: Input the vector data to be converted
:param refer_raster: input reference raster data
:param attribute: Enter the vector field corresponding to the raster value
:param output_raster: output raster data path
:return: None
"""
ds_raster = gdal.Open(refer_raster)
ds_proj = ds_raster.GetProjection() # Projection information
ds_trans = ds_raster.GetGeoTransform() # Affine geographic transformation parameters
ds_width = ds_raster.RasterXSize # Get width/number of columns
ds_height = ds_raster.RasterYSize # Get the height/number of rows
ds_raster = None # Release memory
del ds_raster
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(input_shp, 1)
layer = ds.GetLayer() # Get the layer file object
result = gdal.GetDriverByName('GTiff').Create(output_raster, ds_width, ds_height, bands=1, eType=gdal.GDT_Byte)
result.SetGeoTransform(ds_trans) #Write affine geographical transformation parameters
result.SetProjection(ds_proj) #Write projection information
band = result.GetRasterBand(1)
band.SetNoDataValue(0) # Ignore background value
band.FlushCache() # Clear the data cache
# options = ["ATTRIBUTE=attribute", "CHUNKY SIZE=0", "ALL_TOUCHED=False"]
# Specify the field value in the attribute field of the input vector data to be written into the raster file as a raster value, and the value will be output to all output bands. If this value is specified, the value of the burn_Values parameter can be set to NULL.
# Specify the height of the block in which this operation will be run. The larger the value, the smaller the calculation time required. If this value is not set or is set to 0, the cache size of GDAL is obtained according to the formula: number of bytes occupied by the cache / number of bytes of the scan function. So the value will not exceed the size of the cache.
# Set to TRUE to indicate that all pixels touch the vector center line or polygon, otherwise it is just the center of the polygon or the part selected by the Bresenham algorithm. The default value is FALSE. Simply put, FALSE means eight-direction rasterization, and TRUE means full-path rasterization.
gdal.RasterizeLayer(result, [1], layer, burn_values=[1], options=[f"ATTRIBUTE={attribute}"])
# Vector to raster function. Output raster, band number, layer, raster value (controllable failure), optional parameters (raster value field)
result = None
del result, layer
3. Complete code
There’s not much to say, the comments are very detailed, just call whichever function you use!
# -*- coding: utf-8 -*-
"""
@Time: 2023/11/3 14:37
@Auth: RS lost little book boy
@File: Convert Raster And Vector.py
@IDE:PyCharm
@Purpose: Convert raster data and vector data to each other
"""
import os
from osgeo import gdal, osr, ogr
def Raster_To_Vector(input_raster, output_shp):
"""
:param input_raster: Enter the raster data to be converted
:param output_shp: path of output vector
:return:None
"""
print("Converting raster to polygon...")
# Convert raster to polygon
ds_raster = gdal.Open(input_raster) # Read raster data in the path
band_raster = ds_raster.GetRasterBand(1) # Get the band that needs to be converted to vector
proj_raster = osr.SpatialReference()
proj_raster.ImportFromWkt(ds_raster.GetProjection())
# Assign the projection information of raster data to vector
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(output_shp): # If the file already exists, delete it and do it again
driver.DeleteDataSource(output_shp)
polygon = driver.CreateDataSource(output_shp) # Create data resources
layer_polygon = polygon.CreateLayer("Shp", srs=proj_raster, geom_type=ogr.wkbMultiPolygon) # Create a layer and define multiple faces
new_field = ogr.FieldDefn('value', ogr.OFTReal) # Add a field to the target shp file to store the pixel value of the original raster
layer_polygon.CreateField(new_field)
gdal.FPolygonize(band_raster, None, layer_polygon, 0) # Core function, which performs raster to vector operation
polygon.SyncToDisk()
polygon=None
def Vector_To_Raster(input_shp, refer_raster, attribute, output_raster):
"""
:param input_shp: Input the vector data to be converted
:param refer_raster: input reference raster data
:param attribute: Enter the vector field corresponding to the raster value
:param output_raster: output raster data path
:return: None
"""
ds_raster = gdal.Open(refer_raster)
ds_proj = ds_raster.GetProjection() # Projection information
ds_trans = ds_raster.GetGeoTransform() # Affine geographic transformation parameters
ds_width = ds_raster.RasterXSize # Get width/number of columns
ds_height = ds_raster.RasterYSize # Get the height/number of rows
ds_raster = None # Release memory
del ds_raster
driver = ogr.GetDriverByName("ESRI Shapefile")
ds = driver.Open(input_shp, 1)
layer = ds.GetLayer() # Get the layer file object
result = gdal.GetDriverByName('GTiff').Create(output_raster, ds_width, ds_height, bands=1, eType=gdal.GDT_Byte)
result.SetGeoTransform(ds_trans) #Write affine geographical transformation parameters
result.SetProjection(ds_proj) #Write projection information
band = result.GetRasterBand(1)
band.SetNoDataValue(0) # Ignore background value
band.FlushCache() # Clear the data cache
# options = ["ATTRIBUTE=attribute", "CHUNKY SIZE=0", "ALL_TOUCHED=False"]
# Specify the field value in the attribute field of the input vector data to be written into the raster file as a raster value, and the value will be output to all output bands. If this value is specified, the value of the burn_Values parameter can be set to NULL.
# Specify the height of the block in which this operation will be run. The larger the value, the smaller the calculation time required. If this value is not set or is set to 0, the cache size of GDAL is obtained according to the formula: number of bytes occupied by the cache / number of bytes of the scan function. So the value will not exceed the size of the cache.
# Set to TRUE to indicate that all pixels touch the vector center line or polygon, otherwise it is just the center of the polygon or the part selected by the Bresenham algorithm. The default value is FALSE. Simply put, FALSE means eight-direction rasterization, and TRUE means full-path rasterization.
gdal.RasterizeLayer(result, [1], layer, burn_values=[1], options=[f"ATTRIBUTE={attribute}"])
# Vector to raster function. Output raster, band number, layer, raster value (controllable failure), optional parameters (raster value field)
result = None
del result, layer
if __name__ == "__main__":
shp = r"B:\Personal\Peng Junxi\Deep Learning\CNN\0919 Small Experiment/sample.shp"
refer = r"B:\Personal\Peng Junxi\Deep Learning\CNN\0919 Small Block Experiment\example.tif"
out = r"B:\Personal\Peng Junxi\Deep Learning\CNN\0919 Small Experiment\try.tif"
Attribute = "AREA"
Vector_To_Raster(shp, refer, Attribute, out)
This article mainly shares some codes that I have written while learning Python. Some parts are based on tutorials from predecessors and the official website. If there is any infringement, please contact the author to delete it. If you have any questions, you can leave a message at any time and the blogger will reply in time.