1

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:

enter image description here

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 happens when you try those first two methods? – PolyGeo Jun 25 '16 at 12:10
  • 1
    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. – Sant Jun 25 '16 at 13:03
  • Why can you only use gdal modules, and not console applications? – BradHards Jun 26 '16 at 4:34
  • Because, it's my task. – Sant Jun 26 '16 at 8:57
4

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)? – Sant Jun 26 '16 at 13:50
  • You use a filter with geom to read the bands. – gene Jun 26 '16 at 13:55
  • I understand this) but how use this filter? Have you an example of code? – Sant 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) – Sant Jun 26 '16 at 14:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.