Raster transformation¶
Transform¶
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’
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
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
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
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
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