In this tutorial, we are gonna convert a Polygon and Points vector data into raster format by using GeoPandas and GeoCube in python.
To get full source code for this tutorial, click here. All the related Polygon and Point data used in this tutorial are in repository.
Before rasterizing the vector, we will set up conda environment and install required python GIS packages for this task. Make sure conda is installed on your system. Use following commands to create a conda environment and to install python libraries.
When you install geocube library using conda, it will automatically install related dependencies that includes gdal, shapely, rasterio, geopandas, xarray, rioxarray etc.
And we will also install pygeos library which is used to speed up the vectorized operations in GeoPandas and Shapely.
|
1 2 3 4 |
(base) geoknight@pop-os:~$conda create -n spatial-dev.guru python=3.10 (base) geoknight@pop-os:~$conda activate spatial-dev.guru (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge geocube (spatial-dev.guru) geoknight@pop-os:~$conda install -c conda-forge pygeos |
In this tutorial, we will rasterize both Polygon and Points data.
1. Import Dependencies
Import following libraries in your code. Pandas and Geopandas will used to read csv and shapefile respectively. Geocube library will be used to rasterize the polygon and vector data.
|
1 2 3 4 5 6 |
import pandas as pd import geopandas as gpd from geocube.api.core import make_geocube from geocube.rasterize import rasterize_image from functools import partial from rasterio.enums import MergeAlg |
2. Rasterize Polygons Data
We are using a shapefile for California which is a polygon dataset and stores population value for each feature. The file is in “shps/census_pops.zip” path. We are using geopandas to read this shapefile directly from zip file.
This dataset is in 4326 coordinate system. We will reproject it to 3857 projected coordinate system at line number 7 to make measurement in meters.
After that we will use make_geocube module from geocube library to rasterize the the vector. This module takes geodataframe in vector_data argument. In our case, the geodataframe is the shapefile that we read using geopandas at line number 3. The measurement parameter specifies which column to use from geodataframe to burn the values in pixel. In our case, it is pops column. The resolution parameter specifies the size of of pixel. In our case, it is 500X500 meter per pixel size.
The make_geocube returns the rioxarray object which then can be used to save the raster. At line number 18, we are saving the rasterized data.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
def polygonsToRaster(): # Read California Census Polygon Shapefile cal_census = gpd.read_file('shps/census_pops.zip') cal_census = cal_census.rename(columns={'POP20': 'pops'}) # Reprojecting to 3857 coordinate system cal_census = cal_census.to_crs('EPSG:3857') # Using GeoCube to rasterize the Vector cal_census_raster = make_geocube( vector_data=cal_census, measurements=["pops"], resolution=(-500, 500), fill = 0 ) # Save raster census raster cal_census_raster.rio.to_raster('rasters/cal_census.tiff') |
3. Rasterize Points Data
We are having the uk accidents records for 2005-2007,2009-2011 and 2012-2014 years in csv and compressed format . Each csv file have latitude and longitude coordinates of each accident.
We are reading these csvs using pandas(line number 3,4,5) and merging them into single pandas dataframe(line number 7).
At line number 10 and 11, we changing the column type to float type of Lat and Lon column. After that, we will use geopandas to convert lat long values in Point Geometry and will create a Point type geodataframe using geopandas. See line number 14. This geodataframe is in 4326 coordinate system. We will reproject this geodataframe to 3857 projected coordinate system.
Now we have a point geodataframe of uk_accidents csvs . We will pass this point geodataframe to rasterize it. Rasterization is same as we did for Polygon. The only difference here is that we are passing a rasterize_function parameter in make_geocube. This parameter will let you decide how you want to burn the values in pixels. Here, this parameter, will calculate the number of points intersects with a given pixel of 500×500 meter size and then sum the values given in measurement parameter and finally burn this value in given pixel.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 |
def pointsToRaster(): # Read UK Accidents Point csvs uk_accidents_1 = pd.read_csv('csvs/uk_accidents_2005_to_2007.tar.xz') uk_accidents_2 = pd.read_csv('csvs/uk_accidents_2009_to_2011.tar.xz') uk_accidents_3 = pd.read_csv('csvs/uk_accidents_2012_to_2014.tar.xz') uk_accidents = pd.concat([uk_accidents_1, uk_accidents_2, uk_accidents_3]) # Convert Long Lat into numeric type uk_accidents['Longitude'] = pd.to_numeric(uk_accidents['Longitude']) uk_accidents['Latitude'] = pd.to_numeric(uk_accidents['Latitude']) # Convert Long Lat into Point Geometry uk_accidents = gpd.GeoDataFrame(geometry = gpd.points_from_xy(x=uk_accidents['Longitude'], y=uk_accidents['Latitude'])) uk_accidents = uk_accidents.set_crs('EPSG:4326') # Reprojecting to 3857 coordinate system uk_accidents = uk_accidents.to_crs('EPSG:3857') uk_accidents['value'] = 1 uk_accidents = uk_accidents[uk_accidents.is_valid] uk_accidents = uk_accidents[~uk_accidents.is_empty] # Using GeoCube to rasterize the Vector uk_accidents_raster = make_geocube( vector_data=uk_accidents, measurements=["value"], resolution=(-500, 500), rasterize_function=partial(rasterize_image, merge_alg=MergeAlg.add), fill = 0 ) # save uk_accidents raster uk_accidents_raster.rio.to_raster('rasters/uk_accidents.tiff') |
That’s it. We have successfully rasterized the both Polygon and Point vector data. Geo_Cube library makes it really easy to rasterize the vector data in Python.
Full Source code is below
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 |
import pandas as pd import geopandas as gpd from geocube.api.core import make_geocube from geocube.rasterize import rasterize_image from functools import partial from rasterio.enums import MergeAlg def polygonsToRaster(): # Read California Census Polygon Shapefile cal_census = gpd.read_file('shps/census_pops.zip') cal_census = cal_census.rename(columns={'POP20': 'pops'}) # Reprojecting to 3857 coordinate system cal_census = cal_census.to_crs('EPSG:3857') # Using GeoCube to rasterize the Vector cal_census_raster = make_geocube( vector_data=cal_census, measurements=["pops"], resolution=(-500, 500), fill = 0 ) # Save raster census raster cal_census_raster.rio.to_raster('rasters/cal_census.tiff') def pointsToRaster(): # Read UK Accidents Point csvs uk_accidents_1 = pd.read_csv('csvs/uk_accidents_2005_to_2007.tar.xz') uk_accidents_2 = pd.read_csv('csvs/uk_accidents_2009_to_2011.tar.xz') uk_accidents_3 = pd.read_csv('csvs/uk_accidents_2012_to_2014.tar.xz') uk_accidents = pd.concat([uk_accidents_1, uk_accidents_2, uk_accidents_3]) # Convert Long Lat into numeric type uk_accidents['Longitude'] = pd.to_numeric(uk_accidents['Longitude']) uk_accidents['Latitude'] = pd.to_numeric(uk_accidents['Latitude']) # Convert Long Lat into Point Geometry uk_accidents = gpd.GeoDataFrame(geometry = gpd.points_from_xy(x=uk_accidents['Longitude'], y=uk_accidents['Latitude'])) uk_accidents = uk_accidents.set_crs('EPSG:4326') # Reprojecting to 3857 coordinate system uk_accidents = uk_accidents.to_crs('EPSG:3857') uk_accidents['value'] = 1 uk_accidents = uk_accidents[uk_accidents.is_valid] uk_accidents = uk_accidents[~uk_accidents.is_empty] # Using GeoCube to rasterize the Vector uk_accidents_raster = make_geocube( vector_data=uk_accidents, measurements=["value"], resolution=(-500, 500), rasterize_function=partial(rasterize_image, merge_alg=MergeAlg.add), fill = 0 ) # save uk_accidents raster uk_accidents_raster.rio.to_raster('rasters/uk_accidents.tiff') if __name__ == '__main__': polygonsToRaster() pointsToRaster() |
I hope this tutorial will create a good foundation for you. If you want tutorials on another GIS topic or you have any queries, please send an email at contact@spatial-dev.guru.