3

I have a shapefile that looks like the follows: enter image description here

And here's how the attribute table looks like:

enter image description here

I'm trying to rasterize it with a reference raster and my code is as follows:

input_vector = city_mobile_dir
output_raster = f'../../../data/world_statistics/output_data/city_mobile_gdf_{net}_{iso}.TIF'
ref_raster = f'{source}/world_pop/jpn_ppp_2019_UNadj.tif'
raster_df = gdal.Open(ref_raster, gdal.GA_ReadOnly)
vector_driver = ogr.GetDriverByName('ESRI Shapefile')
vector_df = vector_driver.Open(input_vector, 0)
vector_layer = vector_df.GetLayer()
# fetch number of rows and columns
ncol = raster_df.RasterXSize
nrow = raster_df.RasterYSize

# fetch projection and extent
proj = raster_df.GetProjectionRef()
ext = raster_df.GetGeoTransform()

# create the raster dataset
memory_driver = gdal.GetDriverByName('GTiff')
out_raster_df = memory_driver.Create(output_raster, ncol, nrow, 1, gdal.GDT_Byte)

# set the ROI image's projection and extent to our input raster's projection and extent
out_raster_df.SetProjection(proj)
out_raster_df.SetGeoTransform(ext)

# fill our output band with the 0 blank, no class label, value
band = out_raster_df.GetRasterBand(1)
band.SetNoDataValue(-9999)
band.FlushCache()

# rasterize the shapefile layer to our new dataset
status = gdal.RasterizeLayer(out_raster_df,  
                             [1],  
                             vector_layer,  
                             None, None,  
                             [0],  # burn value 0
                             ['ATTRIBUTE = FID'])

if status != 0:
    print("Failed")
else:
    print("Success")

But no matter what I tried, I kept getting an entire black tiff. I looked through other similar problems on the website but their solutions don't work on mine either. Unfortunately I'm new to GDAL so I might not have been able to identify the mistake I made in my code. Thanks a lot for helping me!

2
  • 1
    FID (feature id) has a special meaning and it is not usually exposed as a normal attribute. I am not sure if it can be used as a burn attribute. And when you burn you use either fixed value or you take the value from an attribute. It seems that you try to do the both. If you want to rasterize, why don't you use gdal.org/programs/gdal_rasterize.html#gdal-rasterize or the python library version of that? See examples about using the latter from the autotests github.com/OSGeo/gdal/blob/master/autotest/utilities/….
    – user30184
    Oct 28 '20 at 7:05
  • 1
    Also if you have a handful of features and fid is between 0 and 20 or so, the image will look almost black on grayscale 0-255 even if you manage to burn fid into pixels.
    – user30184
    Oct 28 '20 at 8:11
1

I think you want something like this:

def create_labels(sample_file, vector_file,
            attribute="FID", otype="MEM", fname_out=""):
    """Rasterise a vector layer `vector_file` using a sample raster file 
    (`sample_file`) to provide extent and projection. You can select an
    attribute in  `vector_file`  (by default, `FID`), and you can either do the
    processing in memory, or dump to disk in any GDAL  supported format.

    Parameters
    ----------
    sample_file : str
        A raster with a suitable geotransform, projection etc.
    vector_file : str
        The vector file you want rasterised
    attribute :str, optional
        The attribute field to rasterise, by default "FID"
    otype : str, optional
        GDAL output type, by default "MEM"  (in memory). Set to `GTIff` for
        file output
    fname_out : str, optional
        The output filename, by default "" (e.g.  in memory raster). 

    Returns
    -------
    ndarray
        An array with the rasterised raster as a numpy array.
    """
    g = gdal.Open(sample_file)
    dst = gdal.GetDriverByName(otype).Create(fname_out,
                                g.RasterXSize, 
                                g.RasterYSize, 1, 
                                gdal.GDT_UInt32)
    dst.SetGeoTransform(g.GetGeoTransform())
    dst.SetSpatialRef(g.GetSpatialRef())
    f  = ogr.Open(vector_file),
    gdal.RasterizeLayer(dst, [1], f.GetLayerByIndex(0),
                        options=[f'ATTRIBUTE={attribute}'])
    if otype != "MEM":
        # Read data and close dataset to flush to disk if
        # not using in-memory array.
        retval = dst.ReadAsArray()
        dst = None
    else:
        return dst.ReadAsArray()

I believe in your code you set the burn value to be 0.

1
  • quite nice code, but I have to say - if you want to add the permanently as Gtiff you have to add ´dst = None´ at the end
    – David
    Feb 15 at 18:00
0

The issue is, that the data of your raster is stored not permanently yet. To do so you have to add out_raster_df = None. With other packages while writing a file you wouldt do this with a with open(...) statement but this is unfortunately not implementet in the GDAL library.

def rasterize(rasterFile, polygonShape, outputPath, burnValue = 1, datatype = gdal.GDT_Int16):
  #input information provided by the raster
  x_res = rasterFile.RasterXSize
  y_res = rasterFile.RasterYSize
  nbands = 1
  
  #safety-net to delete an old file from a previous run
  if os.path.exists(outputPath):
    os.remove(outputPath)
 
  #create an empty raster with the input information
  empty_raster = gdal.GetDriverByName('GTiff').Create(outputPath, x_res, y_res, nbands, datatype, options=['COMPRESS=LZW', 'TILED=YES'])
  empty_raster.SetGeoTransform(rasterFile.GetGeoTransform())
  empty_raster.SetSpatialRef(rasterFile.GetSpatialRef())

  # store the shape information as burn value
  status = gdal.RasterizeLayer(empty_raster, [1], polygonShape.GetLayerByIndex(0), burn_values=[burnValue])
  empty_raster = None # this is the very important line
  return(status)

This function though calls the raster and shape already opened. If you want to call the function with the filepaths simply add rasterFile = gdal.open(rasterpath) and polygonShape = ogr.open(polygonpath).

This function creates a rectangular raster in the size of the input raster. The location of the polygons get the value 1, the rest of the raster has value 0.

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.