From 23bc94a256e2bc26a970e55d2e31bd811ae3e2b8 Mon Sep 17 00:00:00 2001 From: Sean Gillies Date: Mon, 28 Apr 2014 14:39:49 -0600 Subject: [PATCH] Enable NUM_THREADS option for reproject(). Closes #60. Usage: from rasterio.warp import reproject() reproject(..., num_threads=2) In practice, I'm not seeing a speed up for modestly sized images. I think you'd have to get into some very large rasters to see a difference. --- docs/datasets.rst | 25 +++++++++++++++---------- rasterio/_gdal.pxd | 2 ++ rasterio/_warp.pyx | 25 ++++++++++++++++++++++--- rasterio/tests/test_warp.py | 30 ++++++++++++++++++++++++++++++ 4 files changed, 69 insertions(+), 13 deletions(-) diff --git a/docs/datasets.rst b/docs/datasets.rst index ca907de0..d0670e8d 100644 --- a/docs/datasets.rst +++ b/docs/datasets.rst @@ -38,25 +38,16 @@ Attributes In addition to the file-like attributes shown above, a dataset has a number of other read-only attributes that help explain its role in spatial information -systems. +systems. The ``driver`` .. code-block:: pycon >>> dataset.driver u'GTiff' - >>> dataset.shape - (718, 791) >>> dataset.height, dataset.width (718, 791) >>> dataset.shape (718, 791) - >>> dataset.transform - [101985.0, 300.0379266750948, 0.0, 2826915.0, 0.0, -300.041782729805] - >>> dataset.crs - {u'units': u'm', u'no_defs': True, u'ellps': u'WGS84', u'proj': u'utm', u'zone': 18} - -.. code-block:: pycon - >>> dataset.count 3 >>> dataset.dtypes @@ -64,6 +55,20 @@ systems. >>> dataset.nodatavals [0.0, 0.0, 0.0] +What makes geospatial raster datasets different from other raster files is +that their pixels map to regions of the Earth. A dataset has a coordinate +reference system and an affine transformation matrix that maps pixel +coordinates to coordinates in that reference system. + +.. code-block:: pycon + + >>> dataset.crs + {u'units': u'm', u'no_defs': True, u'ellps': u'WGS84', u'proj': u'utm', u'zone': 18} + >>> dataset.transform + [101985.0, 300.0379266750948, 0.0, 2826915.0, 0.0, -300.041782729805] + + + Reading data ------------ diff --git a/rasterio/_gdal.pxd b/rasterio/_gdal.pxd index 41c19364..c7bffc7e 100644 --- a/rasterio/_gdal.pxd +++ b/rasterio/_gdal.pxd @@ -5,11 +5,13 @@ cdef extern from "cpl_conv.h": void * CPLMalloc (size_t) void CPLFree (void *ptr) void CPLSetThreadLocalConfigOption (char *key, char *val) + const char *CPLGetConfigOption (char *, char *) cdef extern from "cpl_string.h": int CSLCount (char **papszStrList) char ** CSLAddNameValue (char **papszStrList, const char *pszName, const char *pszValue) int CSLFindName (char **papszStrList, const char *pszName) + const char * CSLFetchNameValue (char **papszStrList, const char *pszName) char ** CSLSetNameValue (char **list, char *name, char *val) void CSLDestroy (char **list) diff --git a/rasterio/_warp.pyx b/rasterio/_warp.pyx index 392ea4e2..91cfef98 100644 --- a/rasterio/_warp.pyx +++ b/rasterio/_warp.pyx @@ -86,6 +86,10 @@ def reproject( cdef char *dstwkt= NULL cdef const char *proj_c cdef void *osr + cdef char **warp_extras + cdef char *key_c + cdef char *val_c + cdef const char* pszWarpThreads # If the source is an ndarray, we copy to a MEM dataset. # We need a src_transform and src_dst in this case. These will @@ -223,8 +227,23 @@ def reproject( log.debug("Created transformer") psWOptions = _gdal.GDALCreateWarpOptions() - # TODO: get options from kwargs - # psWOptions = GDALCloneWarpOptions( psOptions ); + + warp_extras = psWOptions.papszWarpOptions + for k, v in kwargs.items(): + k, v = k.upper(), str(v).upper() + key_b = k.encode('utf-8') + val_b = v.encode('utf-8') + key_c = key_b + val_c = val_b + warp_extras = _gdal.CSLSetNameValue(warp_extras, key_c, val_c) + + pszWarpThreads = _gdal.CSLFetchNameValue(warp_extras, "NUM_THREADS") + if pszWarpThreads == NULL: + pszWarpThreads = _gdal.CPLGetConfigOption( + "GDAL_NUM_THREADS", "1") + warp_extras = _gdal.CSLSetNameValue( + warp_extras, "NUM_THREADS", pszWarpThreads) + log.debug("Created warp options") psWOptions.eResampleAlg = <_gdal.GDALResampleAlg>resampling @@ -261,7 +280,7 @@ def reproject( log.debug( "Chunk and warp window: %d, %d, %d, %d", 0, 0, cols, rows) - eErr = oWarper.ChunkAndWarpImage(0, 0, cols, rows) + eErr = oWarper.ChunkAndWarpMulti(0, 0, cols, rows) log.debug("Chunked and warped: %d", retval) except Exception, e: diff --git a/rasterio/tests/test_warp.py b/rasterio/tests/test_warp.py index c551ada2..5ced47e5 100644 --- a/rasterio/tests/test_warp.py +++ b/rasterio/tests/test_warp.py @@ -106,3 +106,33 @@ def test_warp_from_to_file(tmpdir): reproject(rasterio.band(src, i), rasterio.band(dst, i)) # subprocess.call(['open', tiffname]) +def test_warp_from_to_file_multi(tmpdir): + """File to file""" + tiffname = str(tmpdir.join('foo.tif')) + with rasterio.open('rasterio/tests/data/RGB.byte.tif') as src: + dst_transform = [-8789636.708, 300.0, 0.0, 2943560.235, 0.0, -300.0] + dst_crs = dict( + proj='merc', + a=6378137, + b=6378137, + lat_ts=0.0, + lon_0=0.0, + x_0=0.0, + y_0=0, + k=1.0, + units='m', + nadgrids='@null', + wktext=True, + no_defs=True) + kwargs = src.meta.copy() + kwargs.update( + transform=dst_transform, + crs=dst_crs) + with rasterio.open(tiffname, 'w', **kwargs) as dst: + for i in (1, 2, 3): + reproject( + rasterio.band(src, i), + rasterio.band(dst, i), + num_threads=2) + # subprocess.call(['open', tiffname]) +