Geographic Information Systems Asked on September 4, 2021
I’m using this GDAL based rasterize function in python in order to create raster from shapefile based on numerical class values.
I am using:
function of GDALMy class attributes are very long numbers :
so I use those clases for the rasterizing as following (the code is from here: https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html):
#Rasterize function
def rasterise_me(raster, vector, attribute,
fname_out="", format="MEM"):
"""Rasterises a vector dataset by attribute to match a given
raster dataset. This functions allows for the raster and vector
to have different projections, and will ensure that the output
is consistent with the input raster.
By default, it returns a handle to an open GDAL dataset that you
can e.g. `ReadAsArray`. If you want to generate a GTiff on disk,
set format to `GTiff` and `fname_out` to a sensible filename.
Parameters
----------
raster: str
The raster filaname used as input. It will not be overwritten.
vector: str
The vector filename
attribute: str
The attribute that you want to rasterize. Ideally, this is
numeric.
fname_out: str, optional
The output filename.
format: str, optional
The output file format, such as GTiff, or whatever else GDAL
understands
"""
# Open input raster file. Need to do this to figure out
# extent, projection & resolution.
gdal.UseExceptions()
g = gdal.Open(raster)
geoT = g.GetGeoTransform()
nx, ny = g.RasterXSize, g.RasterYSize
srs = g.GetProjection()
min_x = min(geoT[0], geoT[0]+nx*geoT[1])
max_x = max(geoT[0], geoT[0]+nx*geoT[1])
min_y = min(geoT[3], geoT[3] + geoT[-1]*ny)
max_y = max(geoT[3], geoT[3] + geoT[-1]*ny)
# Reproject vector to match raster file
vector_tmp = gdal.VectorTranslate("", vector, format="Memory",
dstSRS=srs)
# Do the magic
ds_dst= gdal.Rasterize(fname_out, vector_tmp, attribute=attribute,
outputSRS=srs, xRes=geoT[1], yRes=geoT[-1],
outputBounds=[min_x, min_y, max_x, max_y],
format=format, outputType=gdal.GDT_Int32)
return ds_dst
rasterize=rasterise_me(r"29052019.tif", r"shape.shp", "shape ID",fname_out="raster.tif")
that works and generated a raster, thought the raster seemes to have the same classes but I thought it’s becuasethe class is very big number:
after that I concat this data about the pixels to big table but then when I check the classes in the big table I can see that they are wrong- they are all the same negative number:
I don’t know where do I lose the numbers values or why.
I have tried to check the dtype of this column but it’s int64 before and after the rasterize process.
My end goal: to preserve the original class attributes in the final results
You create the tiff file as Int32 type and obviously your attribute values are too big for Int32. GDAL GeoTIFF driver does not have support for Int64 https://gdal.org/drivers/raster/gtiff.html. You can have a try with Float64. If it does not preserve attribute values accurately you can perhaps modify your attributes. In your sample all values begin with "721…" and it might be possible to cut them out and add back when data are used. A better but more complicated solution would be create a two-band or RGB image and split the Int64 value into separate bands but that would require more Python than just to run gdal.Rasterize.
Correct answer by user30184 on September 4, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP