4

I have been trying to crop out certain area of interest from a raster type .tif image format. I was successfully able to implement that in R, but since shiny doesn't provide much for free, i have to switch from R to python to make a web app. I searched a lot for this problem but couldn't find much. This is the code in R whose functionality in desired:

library(raster)
library(sp)

br<-raster("br.tif")
plot(br)
ps<-as(extent(90.35, 90.5, 25.45,25.55),"SpatialPolygons")
proj4string(ps)<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
aoi_br=crop(br,ps)
plot(aoi_br)

I tried shapely in python, but had trouble installing it, which i then searched and found that it was a problem with many people. I also tried using "pyper" package to call my R script, but that too had troubles.

How can this be implemented?

3 Answers 3

4

shapely is cool but a bit of an overkill if you simply want to crop your raster. rasterio offers the possibility to read 'windows' (doc), which would give you the same results as R's raster::crop().

import rasterio
from affine import Affine

xmin = -88.23
xmax = -88.1
ymin = 20.64
ymax = 20.72

def window_from_extent(xmin, xmax, ymin, ymax, aff):
    col_start, row_start = ~aff * (xmin, ymax)
    col_stop, row_stop = ~aff * (xmax, ymin)
    return ((int(row_start), int(row_stop)), (int(col_start), int(col_stop)))

with rasterio.open('/path/to/raster.tif') as src:
    aff = src.affine
    meta = src.meta.copy()
    window = window_from_extent(xmin, xmax, ymin, ymax, aff)
    # Read croped array
    arr = src.read(1, window=window)
    # Update dataset metadata (if you need it)
    meta.update(height = window[0][1] - window[0][0],
                width = window[1][1] - window[1][0],
                affine = src.window_transform(window))
    meta.pop('transform', None)
6
  • Hey Loic ! Thanks ! It worked very well for me :) I had just a couple of doubts, what would updating dataset metadata do, and what if i don't update it ? Now i used the cropped image array and after processing i have a final array which is to be reprojetced, i used show( arr, transform=src.window_transform(window) ), to see the plot and i am getting exactly the same extent as what i got using my code in R. When i saw the Rasterio documentation, i couldn't find any example that could reproject the matrix data into a .tif file which i could use with the leaflet for representation. Any idea? Commented Jun 21, 2017 at 12:40
  • The array itself has no notion of geolocation. That information is contained in the dictionary accessed with src.meta. After cropping, extent (origin), width, height are different from the original raster, that's why this information has to be updated. This is only useful if you want to perform spatial operation on the array (like a reprojection), want to write it to a new file, or connect it with another spatial dataset (e.g. rasterizing vector file). Commented Jun 21, 2017 at 15:51
  • Yes, exactly, I got the point of updating the meta data. So, does rasterio provide any method to project that processed array into a raster file with correct extent and crs into a .tif file ? Commented Jun 21, 2017 at 15:59
  • (Re)projecting has a very specific meaning in the spatial world. rasterio.warp can do that, though it doesn't sound like that's what you want to do. To simply write the array to a georeferenced .tif, simply run with rasterio.open('new_file.tif', 'w', **meta) as dst: dst.write(arr, 1) Commented Jun 21, 2017 at 16:29
  • Great Loic ! You got right what i actually wanted, the only thing left is the new_file's 'crs': None, while original files meta contained 'crs':CRS({'init': 'epsg:4326'}), even when am passing all the arguments as you mentioned using **meta. Commented Jun 21, 2017 at 16:59
2

1) Shapely is not very difficult to install if you know Python. Otherwise you can use pygeoif ('shapely ultralight') or directly the geojson module.

from shapely.geometry import box
ps = box(90.35,25.45,90.5, 25.55)
print ps.wkt
POLYGON ((90.5 25.45, 90.5 25.55, 90.35 25.55, 90.35 25.45, 90.5 25.45))
print mapping(ps) #GeoJSON format
{'type': 'Polygon', 'coordinates': (((90.5, 25.45), (90.5, 25.55), (90.35, 25.55), (90.35, 25.45), (90.5, 25.45)),)}

2) For cropping, you can use:
- the Python module GDAL/OGR as in Clip a Raster using a Shapefile or Python GDAL/OGR Cookbook: Clip a GeoTiff with Shapefile
- the rasterio module as in the rasterio cookbook

    -
6
  • I tried installing shapely, but it throws an error saying: OSError: [WinError 126] The specified module could not be found Commented Jun 20, 2017 at 18:21
  • How did you install shapely ?
    – gene
    Commented Jun 20, 2017 at 18:41
  • pip install shapely Also when i tried installing rasterio, it throws error saying: error: command 'C:\\Program Files (x86)\\Microsoft Visual Studio 14.0\\VC\\BIN\\cl.exe' failed with exit status 2 Commented Jun 20, 2017 at 18:56
  • You cannot install shapely or rasterio iwith pip if you work on Windows (no C compiler). Use the versions of Christoph Gohlke
    – gene
    Commented Jun 20, 2017 at 19:06
  • Please forgive my poor knowledge in this field, but i have been working and installing packages for python using PyCharm, so when i downloaded a package from Gohlke site, and followed this command after downloading: pip install <file location>\xxxx-win_amd64.whl it says: xxxx-win_amd64.whl is not a supported wheel on this platform. Commented Jun 20, 2017 at 19:22
0

You can use the rasterio.windows.Window() function.

window = rasterio.windows.Window(99, 0, 1000, 490)

with rasterio.open(filepath) as src:
    subset = src.read(1, window=window)

plt.figure(figsize=(16,18.5))
plt.imshow(subset)

In the function you have to insert column-offset, row-offset, width and heigth to get the section you want.

The source can be found here: https://geohackweek.github.io/raster/04-workingwithrasters/

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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