2

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.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Feng Xinmei
  • 29
  • 1
  • 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 https://gdal.org/programs/gdal_rasterize.html#gdal-rasterize or the python library version of that? See examples about using the latter from the autotests https://github.com/OSGeo/gdal/blob/master/autotest/utilities/test_gdal_rasterize_lib.py. – user30184 Oct 28 '20 at 07: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 08:11

2 Answers2

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.

Jose
  • 3,332
  • 20
  • 21
  • 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 '21 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.

David
  • 393
  • 4
  • 12