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

Accept CRS strings.

Browse Source 
Closes #116.
This commit is contained in:
Sean Gillies 2014-08-04 16:40:35 -06:00
parent da7d45b030
commit e8ec9f85fe
4 changed files with 75 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
docs/georeferencing.rst
View File

@ -31,6 +31,14 @@ attribute.
>>> src.crs_wkt
'PROJCS["UTM Zone 18, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32618"]]'
When opening a new file for writing, you may also use a CRS string as an
argument.
.. code-block:: pycon
>>> with rasterio.open('/tmp/foo.tif', 'w', crs='EPSG:3857', ...) as dst:
... # write data to this Web Mercator projection dataset.
Coordinate Transformation
-------------------------

1
rasterio/_gdal.pxd
View File

@ -23,6 +23,7 @@ cdef extern from "ogr_srs_api.h":
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)

41
rasterio/_io.pyx
View File

@ -1453,25 +1453,32 @@ cdef class RasterUpdater(RasterReader):
log.debug("Input CRS: %r", crs)
# EPSG is a special case.
init = crs.get('init')
if init:
auth, val = init.split(':')
if auth.upper() == 'EPSG':
_gdal.OSRImportFromEPSG(osr, int(val))
# Normally, we expect a CRS dict.
if isinstance(crs, dict):
# EPSG is a special case.
init = crs.get('init')
if init:
auth, val = init.split(':')
if auth.upper() == 'EPSG':
_gdal.OSRImportFromEPSG(osr, int(val))
else:
crs['wktext'] = True
for k, v in crs.items():
if v is True or (k in ('no_defs', 'wktext') and v):
params.append("+%s" % k)
else:
params.append("+%s=%s" % (k, v))
proj = " ".join(params)
log.debug("PROJ.4 to be imported: %r", proj)
proj_b = proj.encode('utf-8')
proj_c = proj_b
_gdal.OSRImportFromProj4(osr, proj_c)
# Fall back for CRS strings like "EPSG:3857."
else:
crs['wktext'] = True
for k, v in crs.items():
if v is True or (k in ('no_defs', 'wktext') and v):
params.append("+%s" % k)
else:
params.append("+%s=%s" % (k, v))
proj = " ".join(params)
log.debug("PROJ.4 to be imported: %r", proj)
proj_b = proj.encode('utf-8')
proj_b = crs.encode('utf-8')
proj_c = proj_b
_gdal.OSRImportFromProj4(osr, proj_c)
_gdal.OSRSetFromUserInput(osr, proj_c)
# Fixup, export to WKT, and set the GDAL dataset's projection.
_gdal.OSRFixup(osr)
_gdal.OSRExportToWkt(osr, &wkt)

42
rasterio/tests/test_write.py
View File

@ -107,6 +107,48 @@ def test_write_crs_transform(tmpdir):
transform=transform,
dtype=rasterio.ubyte) as s:
s.write_band(1, a)
assert s.crs == {'init': 'epsg:32618'}
info = subprocess.check_output(["gdalinfo", name]).decode('utf-8')
assert 'PROJCS["UTM Zone 18, Northern Hemisphere",' in info
# make sure that pixel size is nearly the same as transform
# (precision varies slightly by platform)
assert re.search("Pixel Size = \(300.03792\d+,-300.04178\d+\)", info)
def test_write_crs_transform_2(tmpdir):
"""Using 'EPSG:32618' as CRS."""
name = str(tmpdir.join("test_write_crs_transform.tif"))
a = numpy.ones((100, 100), dtype=rasterio.ubyte) * 127
transform = [101985.0, 300.0379266750948, 0.0,
2826915.0, 0.0, -300.041782729805]
with rasterio.open(
name, 'w',
driver='GTiff', width=100, height=100, count=1,
crs='EPSG:32618',
transform=transform,
dtype=rasterio.ubyte) as s:
s.write_band(1, a)
assert s.crs == {'init': 'epsg:32618'}
info = subprocess.check_output(["gdalinfo", name]).decode('utf-8')
assert 'PROJCS["WGS 84 / UTM zone 18N",' in info
# make sure that pixel size is nearly the same as transform
# (precision varies slightly by platform)
assert re.search("Pixel Size = \(300.03792\d+,-300.04178\d+\)", info)
def test_write_crs_transform_3(tmpdir):
"""Using WKT as CRS."""
name = str(tmpdir.join("test_write_crs_transform.tif"))
a = numpy.ones((100, 100), dtype=rasterio.ubyte) * 127
transform = [101985.0, 300.0379266750948, 0.0,
2826915.0, 0.0, -300.041782729805]
crs_wkt = 'PROJCS["UTM Zone 18, Northern Hemisphere",GEOGCS["WGS 84",DATUM["unknown",SPHEROID["WGS84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-75],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]'
with rasterio.open(
name, 'w',
driver='GTiff', width=100, height=100, count=1,
crs=crs_wkt,
transform=transform,
dtype=rasterio.ubyte) as s:
s.write_band(1, a)
assert s.crs == {'init': 'epsg:32618'}
info = subprocess.check_output(["gdalinfo", name]).decode('utf-8')
assert 'PROJCS["UTM Zone 18, Northern Hemisphere",' in info
# make sure that pixel size is nearly the same as transform