-
Notifications
You must be signed in to change notification settings - Fork 522
/
clip.py
138 lines (111 loc) · 4.88 KB
/
clip.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
"""File translation command"""
import logging
import click
from cligj import format_opt
from .helpers import resolve_inout
from . import options
import rasterio
from rasterio.coords import disjoint_bounds
from rasterio.crs import CRS
from rasterio.windows import Window
logger = logging.getLogger(__name__)
# Geographic (default), projected, or Mercator switch.
projection_geographic_opt = click.option(
'--geographic',
'projection',
flag_value='geographic',
help="Bounds in geographic coordinates.")
projection_projected_opt = click.option(
'--projected',
'projection',
flag_value='projected',
default=True,
help="Bounds in input's own projected coordinates (the default).")
# Clip command
@click.command(short_help='Clip a raster to given bounds.')
@click.argument(
'files',
nargs=-1,
type=click.Path(),
required=True,
metavar="INPUT OUTPUT")
@options.output_opt
@options.bounds_opt
@click.option(
'--like',
type=click.Path(exists=True),
help='Raster dataset to use as a template for bounds')
@format_opt
@projection_geographic_opt
@projection_projected_opt
@options.overwrite_opt
@options.creation_options
@click.pass_context
def clip(ctx, files, output, bounds, like, driver, projection,
overwrite, creation_options):
"""Clips a raster using projected or geographic bounds.
\b
$ rio clip input.tif output.tif --bounds xmin ymin xmax ymax
$ rio clip input.tif output.tif --like template.tif
The values of --bounds are presumed to be from the coordinate
reference system of the input dataset unless the --geographic option
is used, in which case the values may be longitude and latitude
bounds. Either JSON, for example "[west, south, east, north]", or
plain text "west south east north" representations of a bounding box
are acceptable.
If using --like, bounds will automatically be transformed to match the
coordinate reference system of the input.
It can also be combined to read bounds of a feature dataset using Fiona:
\b
$ rio clip input.tif output.tif --bounds $(fio info features.shp --bounds)
"""
from rasterio.warp import transform_bounds
with ctx.obj['env']:
output, files = resolve_inout(files=files, output=output, overwrite=overwrite)
input = files[0]
with rasterio.open(input) as src:
if bounds:
if projection == 'geographic':
bounds = transform_bounds(CRS.from_epsg(4326), src.crs, *bounds)
if disjoint_bounds(bounds, src.bounds):
raise click.BadParameter('must overlap the extent of '
'the input raster',
param='--bounds',
param_hint='--bounds')
elif like:
with rasterio.open(like) as template_ds:
bounds = template_ds.bounds
if template_ds.crs != src.crs:
bounds = transform_bounds(template_ds.crs, src.crs,
*bounds)
if disjoint_bounds(bounds, src.bounds):
raise click.BadParameter('must overlap the extent of '
'the input raster',
param='--like',
param_hint='--like')
else:
raise click.UsageError('--bounds or --like required')
bounds_window = src.window(*bounds)
bounds_window = bounds_window.intersection(
Window(0, 0, src.width, src.height))
# Get the window with integer height
# and width that contains the bounds window.
out_window = bounds_window.round_lengths(op='ceil')
height = int(out_window.height)
width = int(out_window.width)
out_kwargs = src.profile
out_kwargs.update({
'driver': driver,
'height': height,
'width': width,
'transform': src.window_transform(out_window)})
out_kwargs.update(**creation_options)
if 'blockxsize' in out_kwargs and out_kwargs['blockxsize'] > width:
del out_kwargs['blockxsize']
logger.warning("Blockxsize removed from creation options to accomodate small output width")
if 'blockysize' in out_kwargs and out_kwargs['blockysize'] > height:
del out_kwargs['blockysize']
logger.warning("Blockysize removed from creation options to accomodate small output height")
with rasterio.open(output, 'w', **out_kwargs) as out:
out.write(src.read(window=out_window,
out_shape=(src.count, height, width)))