Converts the shapefile of a polygon to a numpy array… here is a solution to the problem.
Converts the shapefile of a polygon to a numpy array
I have a polygon shapefile with its geometry, a list of times and a list of values. The data is as follows:
| TIME | AVG | geometry
---- ----+-----+--------------------------------
0 | NaN | NaN | POLYGON((-0.1599375 51.3977, ...
A Nan value indicates that this particular polygon has no value. I want to convert it to a 3D numpy array (based on time) where each cell will be marked with the value of the grid. I’ve tried this, but there are 0s everywhere.
source_ds = ogr. Open("C:/Users/Nathan/Desktop/pl.shp")
source_layer = source_ds. GetLayer()
pixelWidth = pixelHeight = 0.001
x_min, x_max, y_min, y_max = source_layer. GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)
target_ds = gdal. GetDriverByName('GTiff'). Create('temp.tif', cols, rows, 1, gdal. GDT_Byte)
target_ds. SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds. GetRasterBand(1)
NoData_value = 0
band. SetNoDataValue(NoData_value)
band. FlushCache()
gdal. RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=VALUE"])
target_dsSRS = osr. SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds. SetProjection(target_dsSRS.ExportToWkt())
k=gdal. Open('temp.tif'). ReadAsArray()
Know how to do it? Thanks
Solution
I finally understood.
source_ds = ogr. Open("C:/Users/Nathan/Desktop/pl.shp")
source_layer = source_ds. GetLayer()
pixelWidth = pixelHeight = 0.001
x_min, x_max, y_min, y_max = source_layer. GetExtent()
cols = int((x_max - x_min) / pixelHeight)
rows = int((y_max - y_min) / pixelWidth)
target_ds = gdal. GetDriverByName('GTiff'). Create('temp.tif', cols, rows, 1, gdal. GDT_Byte)
target_ds. SetGeoTransform((x_min, pixelWidth, 0, y_min, 0, pixelHeight))
band = target_ds. GetRasterBand(1)
NoData_value = 0
band. SetNoDataValue(NoData_value)
band. FlushCache()
gdal. RasterizeLayer(target_ds, [1], source_layer, options = ["ATTRIBUTE=AVG"])
target_ds = None #this is the line that makes the difference
gdal. Open('temp.tif'). ReadAsArray()