archive-gh-me/rasterio
1
0
Fork 0
mirror of https://github.com/rasterio/rasterio.git synced 2025-12-08 17:36:12 +00:00
Code Issues Projects Releases Wiki Activity

Added crs check to rio rasterize, relaxed bounds error to warning, added OSR utility functions, increased test coverage

Browse Source
This commit is contained in:
Brendan Ward 2015-04-05 22:38:11 -07:00
parent 5bb1c4fe61
commit df72bcf3ee
5 changed files with 159 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
rasterio/_base.pyx
View File

@ -712,3 +712,29 @@ def _transform(src_crs, dst_crs, xs, ys, zs):
_gdal.OSRDestroySpatialReference(src)
_gdal.OSRDestroySpatialReference(dst)
return retval
def is_geographic_crs(crs):
cdef void *osr_crs = _osr_from_crs(crs)
cdef int retval = _gdal.OSRIsGeographic(osr_crs)
_gdal.OSRDestroySpatialReference(osr_crs)
return retval == 1
def is_projected_crs(crs):
cdef void *osr_crs = _osr_from_crs(crs)
cdef int retval = _gdal.OSRIsProjected(osr_crs)
_gdal.OSRDestroySpatialReference(osr_crs)
return retval == 1
def is_same_crs(crs1, crs2):
cdef void *osr_crs1 = _osr_from_crs(crs1)
cdef void *osr_crs2 = _osr_from_crs(crs2)
cdef int retval = _gdal.OSRIsSame(osr_crs1, osr_crs2)
_gdal.OSRDestroySpatialReference(osr_crs1)
_gdal.OSRDestroySpatialReference(osr_crs2)
return retval == 1

19
rasterio/_gdal.pxd
View File

@ -18,23 +18,28 @@ cdef extern from "cpl_string.h":
void CSLDestroy (char **list)
cdef extern from "ogr_srs_api.h":
void * OCTNewCoordinateTransformation (void *source, void *dest)
void OCTDestroyCoordinateTransformation (void *source)
int OCTTransform (void *ct, int nCount, double *x, double *y, double *z)
int OSRAutoIdentifyEPSG (void *srs)
void OSRCleanup ()
void * OSRClone (void *srs)
void OSRDestroySpatialReference (void *srs)
int OSRExportToProj4 (void *srs, char **params)
int OSRExportToWkt (void *srs, char **params)
int OSRImportFromEPSG (void *srs, int code)
int OSRImportFromProj4 (void *srs, char *proj)
int OSRSetFromUserInput (void *srs, char *input)
int OSRAutoIdentifyEPSG (void *srs)
int OSRFixup(void *srs)
const char * OSRGetAuthorityName (void *srs, const char *key)
const char * OSRGetAuthorityCode (void *srs, const char *key)
int OSRImportFromEPSG (void *srs, int code)
int OSRImportFromProj4 (void *srs, char *proj)
int OSRIsGeographic(void *srs)
int OSRIsProjected(void *srs)
int OSRIsSame(void *srs1, void *srs2)
void * OSRNewSpatialReference (char *wkt)
void OSRRelease (void *srs)
void * OCTNewCoordinateTransformation (void *source, void *dest)
void OCTDestroyCoordinateTransformation (void *source)
int OCTTransform (void *ct, int nCount, double *x, double *y, double *z)
int OSRSetFromUserInput (void *srs, char *input)
cdef extern from "gdal.h" nogil:
void GDALAllRegister()

48
rasterio/rio/features.py
View File

@ -165,7 +165,7 @@ def shapes(
'reference system. Pixels assumed to be square if this option is '
'used once, otherwise use: '
'--res pixel_width --res pixel_height')
@click.option('--src_crs', default='EPSG:4326',
@click.option('--src_crs', default=None, #'EPSG:4326'
help='Source coordinate reference system. Limited to EPSG '
'codes for now. Used as output coordinate system if output does '
'not exist or --like option is not used. Default: EPSG:4326')
@ -197,24 +197,22 @@ def rasterize(
If the output raster exists, rio-rasterize will rasterize feature values
into all bands of that raster. The GeoJSON is assumed to be in the same
coordinate reference system as the output.
coordinate reference system as the output unless --src_crs is provided.
--default_value or property values when using --property must be using a
data type valid for the data type of that raster.
If a template raster is provided using the --like option, the affine
transform, coordinate reference system, and data type from that raster will
be used to create the output.
transform and data type from that raster will be used to create the output.
The GeoJSON is assumed to be in the same coordinate reference system unless
--src_crs is provided.
--default_value or property values when using --property must be using a
data type valid for the data type of that raster.
--driver, --bounds, --dimensions, --res, and --src_crs are ignored when
output exists or --like raster is provided
The GeoJSON must be within the bounds of the existing output or --like
raster, or an error will be raised.
--driver, --bounds, --dimensions, and --res are ignored when output exists
or --like raster is provided
If the output does not exist and --like raster is not provided, the input
@ -233,7 +231,7 @@ def rasterize(
added in the future.
"""
import numpy as np
from rasterio._base import is_geographic_crs, is_same_crs
from rasterio.features import rasterize
from rasterio.features import bounds as calculate_bounds
@ -242,6 +240,8 @@ def rasterize(
files = list(files)
output = files.pop()
input = click.open_file(files.pop(0) if files else '-')
has_src_crs = src_crs is not None
src_crs = src_crs or 'EPSG:4326'
# If values are actually meant to be integers, we need to cast them
# as such or rasterize creates floating point outputs
@ -278,10 +278,15 @@ def rasterize(
if os.path.exists(output):
with rasterio.open(output, 'r+') as out:
if disjoint_bounds(geojson_bounds, out.bounds):
raise click.BadParameter('GeoJSON outside bounds of '
if has_src_crs and not is_same_crs(src_crs, out.crs):
raise click.BadParameter('GeoJSON does not match crs of '
'existing output raster',
param=input, param_hint='input')
param='input', param_hint='input')
if disjoint_bounds(geojson_bounds, out.bounds):
warnings.warn('GeoJSON outside bounds of existing output '
'raster. Are they in different coordinate '
'reference systems?')
meta = out.meta.copy()
@ -305,10 +310,16 @@ def rasterize(
else:
if like is not None:
template_ds = rasterio.open(like)
if has_src_crs and not is_same_crs(src_crs, template_ds.crs):
raise click.BadParameter('GeoJSON does not match crs of '
'--like raster',
param='input', param_hint='input')
if disjoint_bounds(geojson_bounds, template_ds.bounds):
raise click.BadParameter('GeoJSON outside bounds of '
'--like raster', param=input,
param_hint='input')
warnings.warn('GeoJSON outside bounds of --like raster. '
'Are they in different coordinate reference '
'systems?')
kwargs = template_ds.meta.copy()
kwargs['count'] = 1
@ -317,11 +328,12 @@ def rasterize(
else:
bounds = bounds or geojson_bounds
if src_crs == 'EPSG:4326':
if is_geographic_crs(src_crs):
if (bounds[0] < -180 or bounds[2] > 180 or
bounds[1] < -80 or bounds[3] > 80):
raise click.BadParameter(
'Bounds are beyond the valid extent for EPSG:4326.',
'Bounds are beyond the valid extent for geographic '
'coordinates.',
ctx, param=bounds, param_hint='--bounds')
if dimensions:

54
tests/test_crs.py
View File

@ -5,6 +5,7 @@ import sys
import rasterio
from rasterio import crs
from rasterio._base import is_geographic_crs, is_projected_crs, is_same_crs
logging.basicConfig(stream=sys.stderr, level=logging.DEBUG)
@ -56,6 +57,15 @@ def test_write_3857(tmpdir):
AUTHORITY["EPSG","3857"]]""" in info.decode('utf-8')
def test_from_epsg():
crs_dict = crs.from_epsg(4326)
assert crs_dict['init'].lower() == 'epsg:4326'
# Test with invalid EPSG code
with pytest.raises(ValueError):
assert crs.from_epsg(0)
def test_bare_parameters():
""" Make sure that bare parameters (e.g., no_defs) are handled properly,
even if they come in with key=True. This covers interaction with pyproj,
@ -64,3 +74,47 @@ def test_bare_parameters():
# Example produced by pyproj
crs_dict = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
assert crs_dict.get('no_defs', False) is True
crs_dict = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=False +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
assert crs_dict.get('no_defs', True) is False
def test_is_geographic():
assert is_geographic_crs({'init': 'EPSG:4326'}) is True
assert is_geographic_crs({'init': 'EPSG:3857'}) is False
wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
assert is_geographic_crs(wgs84_crs) is True
nad27_crs = crs.from_string('+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs')
assert is_geographic_crs(nad27_crs) is True
lcc_crs = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
assert is_geographic_crs(lcc_crs) is False
def test_is_projected():
assert is_projected_crs({'init': 'EPSG:3857'}) is True
assert is_projected_crs({'init': 'EPSG:4326'}) is False
lcc_crs = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
assert is_projected_crs(lcc_crs) is True
wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
assert is_projected_crs(wgs84_crs) is False
def test_is_same_crs():
crs1 = {'init': 'EPSG:4326'}
crs2 = {'init': 'EPSG:3857'}
assert is_same_crs(crs1, crs1) is True
assert is_same_crs(crs1, crs2) is False
wgs84_crs = crs.from_string('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
assert is_same_crs(crs1, wgs84_crs) is True
# Make sure that same projection with different parameter are not equal
lcc_crs1 = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0')
lcc_crs2 = crs.from_string('+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=45 +lat_0=0')
assert is_same_crs(lcc_crs1, lcc_crs2) is False

48
tests/test_rio_features.py
View File

@ -232,6 +232,18 @@ def test_rasterize(tmpdir, runner):
assert (data == 0).sum() == 55
assert (data == 1).sum() == 145
# Test that src_crs is written into new output
output = str(tmpdir.join('test4.tif'))
result = runner.invoke(features.rasterize,
[output,
'--dimensions', 20, 10,
'--src_crs', 'EPSG:3857'
],
input=TEST_MERC_FEATURECOLLECTION)
assert result.exit_code == 0
with rasterio.open(output) as out:
assert out.crs['init'].lower() == 'epsg:3857'
def test_rasterize_existing_output(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@ -266,6 +278,14 @@ def test_rasterize_existing_output(tmpdir, runner):
assert (data == 1).sum() == 125
assert (data == 2).sum() == 20
# Confirm that a different src_crs is rejected, even if a geographic crs
result = runner.invoke(features.rasterize,
[output,
'--res', 0.5,
'--src_crs', 'EPSG:4269'
], input=TEST_FEATURES)
assert result.exit_code == 2
def test_rasterize_like(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@ -286,6 +306,15 @@ def test_rasterize_like(tmpdir, runner):
[output, '--like', 'foo.tif'], input=TEST_FEATURES)
assert result.exit_code == 2
# Test that src_crs different than --like raster crs breaks
output = str(tmpdir.join('test3.tif'))
result = runner.invoke(features.rasterize,
[output,
'--like', 'tests/data/shade.tif',
'--src_crs', 'EPSG:4326'],
input=TEST_FEATURES)
assert result.exit_code == 2
def test_rasterize_property_value(tmpdir, runner):
# Test feature collection property values
@ -319,6 +348,7 @@ def test_rasterize_property_value(tmpdir, runner):
assert (data == 0).sum() == 55
assert (data == 15).sum() == 145
def test_rasterize_out_of_bounds(tmpdir, runner):
output = str(tmpdir.join('test.tif'))
@ -326,18 +356,14 @@ def test_rasterize_out_of_bounds(tmpdir, runner):
result = runner.invoke(features.rasterize,
[output, '--like', 'tests/data/shade.tif'],
input=TEST_FEATURES)
assert result.exit_code == 2
# Test out of bounds of existing output raster (first have to create one)
result = runner.invoke(features.rasterize,
[output,
'--res', 1000,
'--property', 'val',
'--src_crs', 'EPSG:3857'
],
input=TEST_MERC_FEATURECOLLECTION)
assert result.exit_code == 0
assert 'outside bounds' in result.output
assert os.path.exists(output)
with rasterio.open(output) as out:
data = out.read_band(1, masked=False)
assert data.sum() == 0
# Confirm that this does not fail when out of bounds for existing raster
result = runner.invoke(features.rasterize, [output], input=TEST_FEATURES)
assert result.exit_code == 2
assert result.exit_code == 0
assert 'outside bounds' in result.output