2

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?

|improve this question
2

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)
|improve this answer
  • 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? – Akshay Raghav Jun 21 '17 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). – Loïc Dutrieux Jun 21 '17 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 ? – Akshay Raghav Jun 21 '17 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) – Loïc Dutrieux Jun 21 '17 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. – Akshay Raghav Jun 21 '17 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

    -
|improve this answer
  • I tried installing shapely, but it throws an error saying: OSError: [WinError 126] The specified module could not be found – Akshay Raghav Jun 20 '17 at 18:21
  • How did you install shapely ? – gene Jun 20 '17 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 – Akshay Raghav Jun 20 '17 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 Jun 20 '17 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. – Akshay Raghav Jun 20 '17 at 19:22

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.