Clipping GeoTIFF file by GeoJSON?

Question

Input GeoJSON file:

{
  "type": "Feature",
  "properties": {
    "id": 1
  },
  "geometry": {
    "type": "Multipolygon",
    "coordinates": [[[
      [27.397511483867362, 44.161117466010516],
      [27.393924221672666, 44.159751598403503],
      [27.393556666460618, 44.159252063395591],
      [27.393726740035870, 44.158373985750522],
      [27.392040835956994, 44.157378400690988],
      [27.390354358253163, 44.156239034941315],
      [27.390977924658255, 44.152849194060536],
      [27.391438333095618, 44.149298658002031],
      [27.386781918912796, 44.147461728155896],
      [27.384487250437232, 44.146859408403664],
      [27.382636468741264, 44.156671855578281],
      [27.383891699721374, 44.156645049015140],
      [27.384649769913505, 44.157388133683327],
      [27.385547083122507, 44.160232076255667],
      [27.387997850095061, 44.160722084482430],
      [27.390672446485077, 44.161638147279866],
      [27.395361188085396, 44.163429614137918],
      [27.396513835695238, 44.162325787855522],
      [27.397511483867362, 44.161117466010516]
    ]]]
  }
}

Need to clip geotiff file.

Additional terms:

I can use only gdal modules (without using any console utilities).
math, numpy, etc. are welcome
result image must be like this:

UPDATE:

def crop_by_bounding_box(min_x, max_x, min_y, max_y):
xValues = [min_x, max_x]
yValues = [min_y, max_y]
#
# Register Imagine driver and open file
#
driver = gdal.GetDriverByName('GTiff')
driver.Register()
dataset = gdal.Open(filename + '.tif')
if dataset is None:
    print 'Could not open ' + filename + '.tif'
    sys.exit(1)
#
# Getting image dimensions
#
cols = dataset.RasterXSize
rows = dataset.RasterYSize
bands = dataset.RasterCount
#
# Getting georeference info
#
transform = dataset.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = -transform[5]
#
# Computing Point1(i1,j1), Point2(i2,j2)
#
i1 = int((xValues[0] - xOrigin) / pixelWidth)
j1 = int((yOrigin - yValues[0]) / pixelHeight)
i2 = int((xValues[1] - xOrigin) / pixelWidth)
j2 = int((yOrigin - yValues[1]) / pixelHeight)

new_cols = i2 - i1 + 1
new_rows = j2 - j1 + 1
#
# Create  list to store band data in
#
band_list = []
#
# Read in bands and store all the data in bandList
#
for i in range(bands):
    band = dataset.GetRasterBand(i+1) # 1-based index
    data = band.ReadAsArray(i1, j1, new_cols, new_rows)
    band_list.append(data)

new_x = xOrigin + i1 * pixelWidth
new_y = yOrigin - j1 * pixelHeight
new_transform = (new_x, transform[1], transform[2], new_y, transform[4], transform[5])
#
# Create gtif file
#
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(output_file, new_cols, new_rows, 3, gdal.GDT_Byte)
#
# Writting output raster
#
for j in range(bands):
    data = band_list[j]
    dst_ds.GetRasterBand(j+1).WriteArray(data)
#
# Setting extension of output raster
#
dst_ds.SetGeoTransform(new_transform)
wkt = dataset.GetProjection()
#
# Setting spatial reference of output raster
#
srs = osr.SpatialReference()
srs.ImportFromWkt(wkt)
dst_ds.SetProjection( srs.ExportToWkt() )
#
# Close output raster dataset
#
dataset = None
dst_ds = None

def crop_by_geojson(geojson):
    # extract the geometry
    geom = geojson['geometry']
    # convert to json string for ogr
    geom_str = json.dumps(geom)
    # create the ogr geometry
    geom_ogr = ogr.CreateGeometryFromJson(geom_str)
    # export to wkt
    geom_ogr.ExportToWkt()

    #lat_lon_to_utm(lat, lon)
    min_x, max_x, min_y, max_y = bbox(geojson)
    print min_x, max_x, min_y, max_y
    XYmin = lat_lon_to_utm(max_y, min_x)
    XYmax = lat_lon_to_utm(min_y, max_x)
    crop_by_bounding_box(XYmin[0], XYmax[0], XYmin[1], XYmax[1])

Code to explode a GeoJSON geometry's coordinates object and yield coordinate tuples:

def explode(coords):
    for e in coords:
        if isinstance(e, (float, int, long)):
            yield coords
            break
        else:
            for f in explode(e):
                yield f

a, b = zip(*list(explode(geojson['geometry']['coordinates'])))

But anyway I'm stuck with cutting. This way doesn't work. I have an error:

clip = gdalnumeric.choose(mask, (clip, 0)).astype(gdalnumeric.uint8)

ValueError: shape mismatch: objects cannot be broadcast to a single shape

What is the difference with your first part of the question GDAL python cut geotiff image with geojson file ? Try the two methods and make your conclusions. — gene, Jun 25 '16 at 12:23
I just don't know how to convert .geojson file to .shp. If you give me an advice I can try to do first method. About second method: I don't know how to implement clipping by geometry. — Andrii, Jun 25 '16 at 13:03
Why can you only use gdal modules, and not console applications? — BradHards, Jun 26 '16 at 04:34
a GeoJson string is a simple Python dictionary, therefore it is easy to extract the coordinates and iterate to project the points — gene, Jun 28 '16 at 16:42
Ok I already done this. Finally, I got 2 tuple with coordinates in UTM. What I must do next to cut by using this coordinates? — Andrii, Jun 28 '16 at 17:46

score 4 · Answer 1 · edited Apr 13 '17 at 12:34

1)

I just don't know how to convert .geojson file to .shp.

This is a one of the bases of ogr Python. If you have a geometry, it is very easy to convert it to a shapefile

# geojson is GeoJson Polygon
from osgeo import ogr
output = "geojson.shp"
driver = ogr.GetDriverByName('ESRI Shapefile')
if os.path.exists(output):
   driver.DeleteDataSource(output)
ds = driver.CreateDataSource(output)
# you works with polygons
layer = ds.CreateLayer('geojson', geom_type=ogr.wkbPolygon)
feature_def=layer.GetLayerDefn()
feature = ogr.Feature(feature_def)
geom =ogr.CreateGeometryFromJson(geosjon)
feature.SetGeometryDirectly(geom)
layer.CreateFeature(feature)
feature.Destroy()
ds.Destroy()

But why a shapefile if you can do that directly with geom =ogr.CreateGeometryFromJson(geosjon)as in my answer in GDAL python cut geotiff image with geojson file

2)

I don't know how to implement clipping by geometry.

I give you an example in my answer. When you read a raster with data = band.ReadAsArray() the result is a Numpy array with the raster coordinates and not the geographical coordinates.

ds=gdal.Open('test.tif')
cols = ds.RasterXSize
rows = ds.RasterYSize
print rows, cols
741 990 
band=ds.GetRasterBand(1)
array=band.ReadAsArray()
type(array)
numpy.ndarray
array.shape
(741, 990)

Now the challenge is to extract or read a sub array of this result with the limits of a GeoJson geometry or a shapefile and it is a Numpy problem. For that you need to convert first your limits to pixel offset and size (raster coordinates)

After you have two solutions:

read directly the result (clipped band) with a filter (limits of the geometry);
read the entire band and extract the values after.

What is the best ? It is you problem.

New

If I use your example witch is a Feature Object ("type": "Feature"), and not a simple MultiPolygon ("type": "Multipolygon")

import json
from osgeo import ogr
original = {
"type": "Feature",
 "properties": {
"id": 1},
"geometry": {
"type": "Multipolygon",
"coordinates": [[[
  [27.397511483867362, 44.161117466010516],
  ....
  ]]]
  }}
# extract the geometry
geom = original['geometry']
# convert to json string for ogr
geom_str = json.dumps(geom)
# create the ogr geometry
geom_ogr = ogr.CreateGeometryFromJson(geom_str)
# export to wkt
geom_ogr.ExportToWkt()
'MULTIPOLYGON (((27.397511483867362 44.161117466010516,27.393924221672666 44.159751598403503,27.393556666460618 44.159252063395591,27.39372674003587 44.158373985750522,27.392040835956994 44.157378400690988,27.390354358253163 44.156239034941315,27.390977924658255 44.152849194060536,27.391438333095618 44.149298658002031,27.386781918912796 44.147461728155896,27.384487250437232 44.146859408403664,27.382636468741264 44.156671855578281,27.383891699721374 44.15664504901514,27.384649769913505 44.157388133683327,27.385547083122507 44.160232076255667,27.387997850095061 44.16072208448243,27.390672446485077 44.161638147279866,27.395361188085396 44.163429614137918,27.396513835695238 44.162325787855522,27.397511483867362 44.161117466010516)))'

Ok can you explain how I can clip with numpy after geom =ogr.CreateGeometryFromJson(geosjon)? — Andrii, Jun 26 '16 at 13:50
I understand this) but how use this filter? Have you an example of code? — Andrii, Jun 26 '16 at 13:58
in New in my answer in GDAL python cut geotiff image with geojson file — gene, Jun 26 '16 at 14:01
It's for bounding box, not for polygon. And there is no numpy I really need to cut like on last image(polygon) — Andrii, Jun 26 '16 at 14:19
anyway geom =ogr.CreateGeometryFromJson(geosjon) doesn't work
TypeError: in method 'CreateGeometryFromJson', argument 1 of type 'char const *' — Andrii, Jun 26 '16 at 14:51
You have a problem with your version of GDAL or your GeoJson because it works on Mac OS X or Linux Mint with GDAL 1.11.x and GDAL 2.x and in Create Geometry from GeoJSON (Python GDAL/OGR Cookbook ) — gene, Jun 26 '16 at 15:41

Clipping GeoTIFF file by GeoJSON?

1 Answers1

New