forked from perrygeo/python-rasterstats
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathcli.py
More file actions
127 lines (108 loc) · 4.8 KB
/
Copy pathcli.py
File metadata and controls
127 lines (108 loc) · 4.8 KB
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
# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
import click
from rasterstats import zonal_stats, point_query
from rasterstats.io import read_features
import logging
try:
import simplejson as json
except:
import json
SETTINGS = dict(help_option_names=['-h', '--help'])
version = 0.9
@click.command(context_settings=SETTINGS)
@click.argument('input-geojson', type=click.File('r'), default='-')
@click.argument('output-geojson', type=click.File('w'), default='-')
@click.version_option(version=version, message='%(version)s')
@click.option('--raster', '-r', required=True, type=click.Path(exists=True))
@click.option('--all-touched/--no-all-touched', default=False)
@click.option('--band', type=int, default=1)
@click.option('--categorical/--no-categorical', default=False)
@click.option('--indent', type=int, default=None)
@click.option('--info/--no-info', default=False)
@click.option('--nodata', type=int, default=None)
@click.option('--prefix', type=str, default='_')
@click.option('--stats', type=str, default=None)
def zonalstats(input_geojson, raster, output_geojson, all_touched, band, categorical,
indent, info, nodata, prefix, stats):
'''zonalstats generates summary statistics of geospatial raster datasets
based on vector features.
The input and output arguments of zonalstats should be valid GeoJSON FeatureCollections. The output GeoJSON will be mostly unchanged but have additional properties per feature describing the summary statistics (min, max, mean, etc.) of the underlying raster dataset. The input and output arguments default to stdin and stdout but can also be file paths.
The raster is specified by the required -r/--raster argument.
Example, calculate rainfall stats for each state and output to file:
\b
zonalstats states.geojson -r rainfall.tif > mean_rainfall_by_state.geojson
'''
if info:
logging.basicConfig(level=logging.INFO)
mapping = json.loads(input_geojson.read())
input_geojson.close()
try:
if mapping['type'] == "FeatureCollection":
feature_collection = mapping
else:
feature_collection = {'type': 'FeatureCollection'}
features = read_features(mapping)
except (AssertionError, KeyError):
raise ValueError("input_geojson must be valid GeoJSON")
if stats is not None:
stats = stats.split(" ")
if 'all' in [x.lower() for x in stats]:
stats = "ALL"
zonal_results = zonal_stats(
features,
raster,
all_touched=all_touched,
band_num=band,
categorical=categorical,
nodata=nodata,
stats=stats,
prefix=prefix,
geojson_out=True)
feature_collection['features'] = zonal_results
output_geojson.write(json.dumps(feature_collection, indent=indent))
output_geojson.write("\n")
@click.command(context_settings=SETTINGS)
@click.argument('input-geojson', type=click.File('r'), default='-')
@click.argument('output-geojson', type=click.File('w'), default='-')
@click.version_option(version=version, message='%(version)s')
@click.option('--raster', '-r', required=True, type=click.Path(exists=True))
@click.option('--band', type=int, default=1)
@click.option('--nodata', type=int, default=None)
@click.option('--indent', type=int, default=None)
@click.option('--interpolate', type=str, default='bilinear')
@click.option('--property-name', type=str, default='value')
def pointquery(input_geojson, raster, output_geojson, band, indent, nodata,
interpolate, property_name):
"""
Queries the raster values at the points of the input GeoJSON Features.
The raster values are added to the features properties and output as GeoJSON
Feature Collection.
If the Features are Points, the point geometery is used.
For other Feauture types, all of the verticies of the geometry will be queried.
For example, you can provide a linestring and get the profile along the line
if the verticies are spaced properly.
You can use either bilinear (default) or nearest neighbor interpolation.
"""
mapping = json.loads(input_geojson.read())
input_geojson.close()
try:
if mapping['type'] == "FeatureCollection":
feature_collection = mapping
else:
feature_collection = {'type': 'FeatureCollection'}
features = read_features(mapping)
except (AssertionError, KeyError):
raise ValueError("input_geojson must be valid GeoJSON")
results = point_query(
features,
raster,
band=band,
nodata=nodata,
interpolate=interpolate,
property_name=property_name,
geojson_out=True)
feature_collection['features'] = results
output_geojson.write(json.dumps(feature_collection, indent=indent))
output_geojson.write("\n")