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
Improve this question
5
  • please edit your question with further questions, and don't post follow-ups as answers. Commented Dec 2, 2020 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) Commented Dec 2, 2020 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. Commented Dec 4, 2020 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! Commented Dec 7, 2020 at 17:07
  • It is easy to do that with just GDAL. See here. Commented Dec 15, 2020 at 19:04

1 Answer 1

Reset to default
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.

Improve this answer
4
  • 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.. Commented Dec 8, 2020 at 5:26
  • oops, sorry, I made a mistake, fixed it. Commented Dec 8, 2020 at 14:54
  • Hii, @Gevaert, Thanks for the suggestions. I was able to generate the raster from vector_ESRI_Shapefile. Commented Dec 9, 2020 at 6: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? Commented Dec 9, 2020 at 6:20

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.