In [63]:

    

%matplotlib inline
import matplotlib.pyplot as plt
import skimage.draw
import numpy as np
import shapely.geometry
import osgeo.ogr
import osgeo.gdal
import shapely.wkb


    

In [64]:

    

# example:
# convert from numpy arrays
y = np.array([4000,4000, 4500, 5000,5000])
x = np.array([1000, 2000, 1500, 2000, 1000])
# to shapely
poly = shapely.geometry.Polygon(
    shell=shapely.geometry.LinearRing(
        coordinates=np.c_[x,y]
    )
)
# to wkb
wkb = poly.wkb
# to ogr
geom = osgeo.ogr.CreateGeometryFromWkb(wkb)
# to json
geom_json = geom.ExportToJson()
# and back to ogr
geom = osgeo.ogr.CreateGeometryFromJson(geom_json)
# and back to wkb
wkb = geom.ExportToWkb()
# and back to shapely
poly = shapely.wkb.loads(wkb)
# and back to numpy vectors
x, y = np.array(poly.exterior.coords).T
x, y


    

    
Out[64]:





(array([ 1000.,  2000.,  1500.,  2000.,  1000.,  1000.]),
 array([ 4000.,  4000.,  4500.,  5000.,  5000.,  4000.]))

In [65]:

    

from shapely.geometry import mapping, shape
import json
geom_json = json.dumps(mapping(poly))
shape(json.loads(geom_json))


    

    
Out[65]:





<shapely.geometry.polygon.Polygon at 0x10e994490>

In [66]:

    

filename = '/Users/baart_f/models/duifp/subgrid/crop_Duifp.tif'
#driver = osgeo.gdal.GetDriverByName('Gtiff')
dataset = osgeo.gdal.Open(filename)
band = dataset.GetRasterBand(1)
x0, dxx, dxy, y0, dyx, dyy = dataset.GetGeoTransform()
xrange = (x0, x0 + dxx * dataset.RasterXSize)
yrange = (y0, y0 + dyy * dataset.RasterYSize)
arr = band.ReadAsArray()
arr = np.ma.masked_equal(arr, band.GetNoDataValue())


    

In [71]:

    

def draw_shape_on_raster(geojson, raster, value):
    """
    draw the polygon geojson on the raster with value=value, inline
    """
    
    geom = shape(json.loads(geojson))
    # and back to numpy vectors
    if geom.type == 'Polygon':
        coords = np.array(geom.exterior.coords)
        x, y = coords.T
        rr, cc = skimage.draw.polygon(x=x, y=y, shape=raster.shape)
    else:
        raise ValueError("Can only draw polygon exteriors")
    raster[rr, cc] = value
draw_shape_on_raster(geom_json, arr, -2)


    

In [72]:

    

plt.matshow(arr, cmap='Set2')
plt.plot(x, y, 'k-', alpha=0.3)


    

    
Out[72]:





[<matplotlib.lines.Line2D at 0x10e98e6d0>]






    

In [70]:

    

y = np.array([1, 1, 3, 3])
x = np.array([1, 3, 3, 1])
# to shapely
poly = shapely.geometry.Polygon(
    shell=shapely.geometry.LinearRing(
        coordinates=np.c_[x,y]
    )
)
# to wkb
wkb = poly.wkb
# to ogr
geom = osgeo.ogr.CreateGeometryFromWkb(wkb)
# to json
geom_json = geom.ExportToJson()
raster = np.zeros((4, 4), dtype='int')
draw_shape_on_raster(geom_json, raster, 1)
raster.flatten()


    

    
Out[70]:





array([0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0])

In [51]:

    

import numpy.testing as npt


    

In [53]:

    

npt.assert_array_equal?


    

In [54]:

    

geom_json


    

    
Out[54]:





'{ "type": "Polygon", "coordinates": [ [ [ 1.0, 1.0 ], [ 3.0, 1.0 ], [ 3.0, 3.0 ], [ 1.0, 3.0 ], [ 1.0, 1.0 ] ] ] }'

In [ ]: