This is part of a training course on intermediate Geospatial Programming with Python with a focus on Raster and Earth Observation Data. If you would like to know more please visit www.acgeospatial.co.uk/training or email info@acgeospatial.co.uk

In [2]:
#!/usr/bin/python
# -*- coding: utf-8 -*-
## Written by Andrew Cutts
## I am indebted to code from the GDAL OGR Cookbook
## https://pcjericks.github.io/py-gdalogr-cookbook/


## first up lets import all the python libraries we need
from osgeo import gdal,osr,ogr
import os

We will build two helper functions, firstly to build the boundary of a raster file. Typically this would be a satellite image but it does not have to be

In [3]:
def build_bounds(coords, srs, name):
    ## Function to build geojson polygon from ordered coordinate set, projection and output name
    ####
    print "creating poly"
    # Create ring
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for i in range (0,len(coords)):
        ring.AddPoint(coords[i][0], coords[i][1])

    # Create polygon
    poly = ogr.Geometry(ogr.wkbPolygon)
    poly.AddGeometry(ring)


    # Create the output Driver
    outDriver = ogr.GetDriverByName('GeoJSON')

    # Create the output GeoJSON
    outDataSource = outDriver.CreateDataSource(name)
    outLayer = outDataSource.CreateLayer(name, srs, geom_type=ogr.wkbPolygon )

    # Get the output Layer's Feature Definition
    featureDefn = outLayer.GetLayerDefn()

    # create a new feature
    outFeature = ogr.Feature(featureDefn)

    # Set new geometry
    outFeature.SetGeometry(poly)

    # Add new feature to output Layer
    outLayer.CreateFeature(outFeature)

    # dereference the feature
    outFeature = None

    # Save and close DataSources
    outDataSource = None
    
    return None

Secondly we will read in our Raster dataset and extract the information about its projection & extents

In [4]:
def read_raster(raster_in):	
    ## Open raster_in as ds (abbreviation for dataset)
    ds=gdal.Open(raster_in)
    
    ## Get the projection of the raster and print out its EPSG code
    prj=ds.GetProjection()
    srs=osr.SpatialReference(wkt=prj)
    print "EPSG code: ", srs.GetAttrValue('authority', 1)

    ## using GetGeoTransform we can get the upper left X and upper left y coordinates
    ulx, xres, xskew, uly, yskew, yres  = ds.GetGeoTransform()
    
    ## Calculate lower right x and lower right y we have the coordinates to build a polygon
    ## these values will be returned and be inputs into build bounds
    lrx = ulx + (ds.RasterXSize * xres)
    lry = uly + (ds.RasterYSize * yres)

    ## rows and columns of the imagery (if needed)
    cols = ds.RasterXSize
    rows = ds.RasterYSize

    ## print it all out to command line
    print "Number of columns: " + str(cols)
    print "Number of rows: " + str(rows)

    print "___"
    print "upper left x, upper left y, lower right x, lower right y:"
    print ulx
    print uly
    print lrx
    print lry

    coords = [(ulx,lry), (ulx,uly), (lrx,uly), (lrx,lry)]

    print "coords"
    print coords
    return coords, srs
In [5]:
### lets call the functions
image = ".../output_images/L2A_T30UXB_20170102T111442_TCI_60m.jp2" ## this is a level2a Sentinel 2a image @60m, true colour RGB
coords, srs = read_raster(image)
EPSG code:  32630
Number of columns: 1830
Number of rows: 1830
___
upper left x, upper left y, lower right x, lower right y:
600000.0
5700000.0
709800.0
5590200.0
coords
[(600000.0, 5590200.0), (600000.0, 5700000.0), (709800.0, 5700000.0), (709800.0, 5590200.0)]
In [6]:
## lets now build a geojson file for the footprint

## specify the directory to write
os.chdir('.../output_images') ## specifiy your output directory

## specifiy the name of the geojson file to write to
name = 'L2A_T30UXB_20170102T111442_TCI_60m_boundry.geojson'

## call the first function passing in coords, srs (projection) and the output name
build_bounds(coords, srs, name)
creating poly
In [ ]: