.. _trastertransformation:
    
Raster transformation
_____________________


Transform
+++++++++

.. code-block:: python
    :emphasize-lines: 4

    from girs.rast.raster import RasterReader
    r0 = RasterReader('D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif')
    print r0.get_coordinate_system()
    r1 = r0.transform(epsg=4326)
    print r1.get_coordinate_system()

    
Output::

    PROJCS["WGS 84 / UTM zone 23S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID[...
    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,]

    
Three types of transform can be given in `transform`: 'proj4', 'wkt', and 'epsg'

.. code-block:: python
    :emphasize-lines: 10, 13, 16

    from girs.srs import srs_from_wkt, export
    from girs.rast.proc import composite
    rasters = [
        'D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif',
        'D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.11.tif']
    r0 = composite(*rasters)
    srs0 = r0.get_coordinate_system()
    srs = srs_from_wkt(srs0)

    r1 = r0.transform(epsg=4326)
    srs1 = r1.get_coordinate_system()

    r2 = r1.transform(proj4=export(srs, 'proj4'))
    srs2 = r2.get_coordinate_system()

    r3 = r1.transform(wkt=export(srs, 'wkt'))
    srs3 = r3.get_coordinate_system()

    print srs0
    print srs1
    print srs2
    print srs3

Output::

    PROJCS["WGS 84 / UTM zone 23S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID...
    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,...
    PROJCS["UTM Zone 23, Southern Hemisphere",GEOGCS["WGS 84",DATUM["WGS_1984",...
    PROJCS["WGS 84 / UTM zone 23S",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID[...]



Pixel and world coordinates
+++++++++++++++++++++++++++


**Pixel to world coordinates**

.. code-block:: python
    :emphasize-lines: 5, 6

    from girs.rast.raster import RasterReader
    r = RasterReader('D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif')
    dx, dy = r.get_pixel_size()
    n_cols, n_rows = r.get_raster_size()
    x0, y0 = r.pixel_to_world(0, 0)
    xn, yn = r.pixel_to_world(n_cols, n_rows)
    env = r.get_extent()
    print 'min x: {:2.8}, max x: {:2.8}, min y: {:2.8}, max y: {:2.8}'.format(*env)
    print
    print 'min x: {:2.8}  {:2.8}'.format(x0, env[0])
    print 'max x: {:2.8}  {:2.8}'.format(xn, env[1])
    print 'min y: {:2.8}  {:2.8}'.format(yn, env[2])
    print 'max y: {:2.8}  {:2.8}'.format(y0, env[3])
    print
    print 'max x: {:2.8}  {:2.8}'.format(xn, x0 + dx * n_cols)
    print 'min y: {:2.8}  {:2.8}'.format(yn, y0 - dy * n_rows)


Output::

    min x: 252029.43, max x: 718719.28, min y: 7632902.4, max y: 8017235.2

    min x: 252029.43  252029.43
    max x: 718719.28  718719.28
    min y: 7632902.4  7632902.4
    max y: 8017235.2  8017235.2

    max x: 718719.28  718719.28
    min y: 7632902.4  7632902.4
    
**World to pixel coordinates**

.. code-block:: python
    :emphasize-lines: 10

    import sys
    from girs.rast.raster import RasterReader, world_to_pixel
    r = RasterReader('D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif')
    dx, dy = r.get_pixel_size()
    x0, y0 = r.pixel_to_world(0, 0)
    print 'sys.float_info.epsilon: {}'.format(sys.float_info.epsilon)
    x1 = x0 * (1 + sys.float_info.epsilon)
    y1 = y0 * (1 - sys.float_info.epsilon)
    for i in range(max(r.get_raster_size())):
        print r.world_to_pixel(x0+i*dx, y0-i*dy), r.world_to_pixel(x1+i*dx, y1-i*dy)

Output::

    sys.float_info.epsilon: 2.22044604925e-16
    (0, 0) (0, 0)
    (0, 1) (1, 1)
    (1, 1) (2, 2)
    (3, 3) (3, 3)
    (4, 3) (4, 4)
    (5, 4) (5, 5)
    (6, 6) (6, 6)
    (7, 6) (7, 7)
    (8, 8) (8, 8)
    (9, 9) (9, 9)
    (10, 9) (10, 10)
    (11, 11) (11, 11)
    (12, 11) (12, 12)
    (13, 13) (13, 13)
    (14, 14) (14, 14)
    (15, 14) (15, 15)
    (16, 16) (16, 16)

    
Extent
++++++

**Extent in world coordinates**

.. code-block:: python
    :emphasize-lines: 5

    from girs.rast.raster import RasterReader, extent_world_to_pixel
    r = RasterReader('D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif')
    xmin, xmax, ymin, ymax = r.get_extent(scale=-0.2)
    gt0 = r.get_geotransform()
    rs1, gt1 = extent_world_to_pixel(gt0, xmin, xmax, ymin, ymax)
    print r.get_raster_size()
    print rs1
    print gt0
    print gt1

Output::

    (17, 14)
    (3, 13, 2, 11)
    (252029.43217803087, 27452.34389348292, 0.0, 8017235.206348279, 0.0, -27452.34389348292)
    [334386.46385847963, 27452.34389348292, 0.0, 7962330.5185613129, 0.0, -27452.34389348292]

**Extent union and intersection of rasters in world coordinates**

.. code-block:: python
    :emphasize-lines: 12, 14

    from girs.rast.raster import RasterReader, RasterWriter, rasters_get_extent
    r0 = RasterReader('D:/tmp/girs/chirps_basin/chirps-v2.0.2016.01.10.tif')
    x_min, y_max = r0.pixel_to_world(0, 0)
    gt = list(r0.get_geotransform())
    gt[0], gt[3] = x_min + 100000, y_max - 100000
    p = r0.get_parameters()
    p.geo_trans = gt
    r1 = RasterWriter(p)
    print 'Initial raster extent:   {:.8}, {:.8}, {:.8}, {:.8}'.format(*r0.get_extent())
    print 'Initial raster extent:   {:.8}, {:.8}, {:.8}, {:.8}'.format(*r1.get_extent())
    print 'Intersection of extents: {:.8}, {:.8}, {:.8}, {:.8}'.format(
        *rasters_get_extent([r0, r1], extent_type='intersection'))
    print 'Union of extents:        {:.8}, {:.8}, {:.8}, {:.8}'.format(
        *rasters_get_extent([r0, r1], extent_type='union'))


Output::

    Initial raster extent:   252029.43, 718719.28, 7632902.4, 8017235.2
    Initial raster extent:   352029.43, 818719.28, 7532902.4, 7917235.2
    Intersection of extents: 352029.43, 718719.28, 7632902.4, 7917235.2
    Union of extents:        252029.43, 818719.28, 7532902.4, 8017235.2