0

Followed the steps from this link to create raster from shapefile using GDAL

Convert an OGR File to a Raster

but it is returning an empty raster. Can someone tell how to define a particular column from the shapefile data to burn?

The code is as follows

import gdal
from osgeo import osr
from osgeo import ogr

raster_path = 'D:\SWAT_HRU_Test\SWAT_Bihar\Bihar_shapefile_python\Full_HRU_GCS_output.tif'


shapefile = 'D:\SWAT_HRU_Test\SWAT_Bihar\Bihar_shapefile_python\Full_HRU_GCS.shp'
source_ds = ogr.Open(shapefile)
source_layer = source_ds.GetLayer()
pixelWidth = pixelHeight = 0.01 
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(raster_path, cols, rows, 1, gdal.GDT_Float32) 
target_ds.SetGeoTransform((x_min, pixelWidth, 0, y_max, 0, -pixelHeight))
target_dsSRS = osr.SpatialReference()
target_dsSRS.ImportFromEPSG(4326)
target_ds.SetProjection(target_dsSRS.ExportToWkt())
band = target_ds.GetRasterBand(1) 
band.SetNoDataValue(-9999) 
gdal.RasterizeLayer(target_ds, [1], source_layer, options=["ATTRIBUTE = HRUGIS"])


this HRUGIS is a column in the shapefile on the basis of which i want to create the raster but it is returning a raster with all values 255


target_ds = None
Vince
  • 20,017
  • 15
  • 45
  • 64
Abhilash Singh Chauhan
  • 662
  • 4
  • 19
  • please edit your question with further questions, and don't post follow-ups as answers. – Fee Dec 02 '20 at 16:43
  • Geocube has quite some dependencies, that may not immediately work. If you run conda create GeoCube --strict-channel-priority -c conda-forge geocube <others, like gdal> and activate the environment, it should work (since Geocube depends on rasterio, xarray and geopandas, those will be installed as well) – Fee Dec 02 '20 at 16:45
  • @GevaertJoep, Respected Sir, since I'm new to stack overflow, I don't know how to 'edit question with further questions'. I wanted to attach the image of code, where I got stuck with some error. So have to post them as answers. – Abhilash Singh Chauhan Dec 04 '20 at 20:59
  • you can edit your question and add anything you want. There should be an 'edit' button somewhere. On another note: Is jupyter notebook installed in the new environment and did you start jupyter from the new conda environment? so: conda activate <env-name>; jupyter-notebook alternatively, follow this guide for jupyter and conda. good luck! – Fee Dec 07 '20 at 17:07
  • It is easy to do that with just GDAL. See here. – Jose Dec 15 '20 at 19:04

1 Answers1

0

edit: Since GeoCube has quite some (nontrivial) dependencies, I recommend creating a new environment with: conda create -n <your_name> -c conda-forge geocube <whatever_else>, optionally also --strict-channel-priority

There is a great module: GeoCube, that does exactly what you want. It is somewhat more high-level, but uses GDAL internally I believe. Then, it becomes as easy as:

from geocube.api.core import make_geocube
import rioxarray
import geopandas as gpd

raster_path = 'your path'
vector_path = 'other path'


source_ds = gpd.read_file(vector_path)
DEM = rioxarray.open_rasterio(raster_path,masked=True)
Cube = make_geocube(vector_data=source_ds['columns1','columns2'],like=DEM)


Cube.columns1.plot()

alternatively, you can leave out the ['columns1','columns2'] to rasterize all attributes of the shapefile

sidenote: shapefile is the ESRI standard for storing vector formats, so I prefer to call it vector inside a program/script, since you can derive them from anywhere.

Fee
  • 1,158
  • 4
  • 14
  • Hii @Gevaert, can you please tell, that in the line Cube = make_geocube(vector_data['columns1','columns2']=LUMap,like=DEM), what is the variable =LUMap, sounds like LandUse map, if so, how we have to define this above? do we have to give its path in the template raster above ? but i think, you have used the DEM raster as template in the code. Thanks.. – Abhilash Singh Chauhan Dec 08 '20 at 05:26
  • oops, sorry, I made a mistake, fixed it. – Fee Dec 08 '20 at 14:54
  • Hii, @Gevaert, Thanks for the suggestions. I was able to generate the raster from vector_ESRI_Shapefile. – Abhilash Singh Chauhan Dec 09 '20 at 06:20
  • But while comparing it with raster generated via ArcGIS using 'feature to raster" tool. Like for making soil moisture raster from shapefile using ArcGIS, if the value in HRU polygon is 185 mm, raster by ArcGIS is showing values in near range, and for large HRU in geocube raster, values were in range of 175-195 too, not much deviation), but the smaller HRUs features were lost, like if there is some small HRU with a value of soil moisture 45 mm, then the geocube raster is overestimating it, showing the value near 110mm. Why so? – Abhilash Singh Chauhan Dec 09 '20 at 06:20