gdal-py-snippets

This repository presents a wide range of different-purpose small examples using the gdal-python bindings (both gdal and ogr). It uses a fair amount of numpy as well, along with other third party packages such as pandas.

The programs presented in this repository were developed using Python 3.6.5 and the following packages versions:

• gdal version: 2.3.3
• numpy version: 1.14.3
• pandas version: 0.24.2

Data

All the snippets in this repository use open data. Here is a table with the description of each dataset used^* and its sources:

Dataset	Description	Type	Extension	Source
COL_msk_alt	Elevation for Colombia (country mask)	Raster	.vrt	DIVA_GIS
wc2.0_10m_prec_01	Average January global precipitation from 1970 to 2000.	Raster	.tif	WorldClim
COL_rails	Colombian railroads	Vector	.shp	DIVA_GIS
ne_110m_admin_0_countries	Administrative boundaries (countries) of the world.	Vector	.shp	Natural Earth

^{*I do not own any of the datasets here presented.}

Snippets

Here is a brief explanation for each snippet:

1. attribute_table_to_csv

Saves a shapefile's attribute table to a given csv file using ogr and pandas. The program iterates through each feature in the shapefile and then iterates through each one of its fields, storing the values in a DataFrame. The DataFrame is then converted to a csv file.

2. compute_proximity

Computes proximity (euclidean distance) in pixels from each cell to a target value using the gdal.ComputeProximity() function. The input raster is the output from rasterize.

3. euclidean_distance

Computes the euclidean distance (in pixels) from each cell to a target value. Similar to compute_proximity but instead of using gdal's built-in function an own implementation of the distance computation is presented. In this case, the target values are the rasterized lines produced in rasterize. This snippet uses numpy to vectorize the distance computation (following Pythagoras' theorem). It creates a distance matrix to each target value (i.e. each value different from NoData). Then, the minimum distance value for each cell is retrieved and a single distance matrix is returned.

4. extract_lowest_cell_values

Extracts an arbitrary percentage (5% in this case) of the lowest cell values in the raster.

5. get_fields_info

Loops through each field in a shapefile and prints information about them.

6. get_pixel_coordinates

Creates two grids with the corresponding x and y coordinates of each pixel center in a raster.

7. get_pixel_values

Gets an array of pixel values correspondent to a set of points within a raster's extent. Points are presented as a tuple of 1D arrays with x and y coordinates. It converts coordinates to array indices and retrieve the values by indexing the data array.

7. raster_up_sample^*

Resamples (up-samples) a raster to a higher resolution. In order to do this it reads the original GeoTransform and changes the pixel width and pixel height to be proportionally smaller. Then it reads the data from the original dataset specifying buf_xsize and buf_ysize so the array where the data is going to be stored fits the new dimensions. When specifying a bigger buffer than the original raster dimensions in the ReadAsArray() method, original values are repeated to fit the new dimensions.

^{*This snippet was slightly adapted from one example presented by Chris Garrard in Geoprocessing with Python's 9th chapter.}

8. rasterize

Rasterizes a shapefile. It uses ogr to open a shapefile and get its extent. Then it creates an empty raster with the same extent as the shapefile and an arbitrary value for the pixel resolution. It uses the gdal.RasterizeLayer() function to burn the raster´s band (i.e. first and only band) with a 1 where the shapefile presents a feature. The rest is set as NoData.