Geografic-Information-System/scripts/r.tileset/r.tileset.py

#!/usr/bin/env python

############################################################################
#
# MODULE:       r.tileset
#
# AUTHOR(S):    Cedric Shock
#               Updated for GRASS7 by Martin Landa, 2009
#
# PURPOSE:      To produce tilings of regions in other projections.
#
# COPYRIGHT:    (C) 2006-2009 by Cedric Shoc, Martin Landa, and GRASS development team
#
#               This program is free software under the GNU General
#               Public License (>=v2). Read the file COPYING that
#               comes with GRASS for details.
#
#############################################################################

# Bugs:
#  Does not know about meridians in projections. However, unlike the usual
#   hack used to find meridians, this code is perfectly happy with arbitrary
#   rotations and flips
#  The following are planned to be fixed in a future version, if it turns out
#   to be necessary for someone:
#   Does not generate optimal tilings. This means that between an appropriate
#    projection for a region and latitude longitude projection, in the
#    degenerate case, it may create tiles demanding up to twice the necessary
#    information. Requesting data from cylindrical projections near their poles
#    results in divergence. You really don't want to use source data from
#    someone who put it in a cylindrical projection near a pole, do you?
#   Not generating "optimal" tilings has another side effect; the sanity
#    of the destination region will not carry over to generating tiles of
#    realistic aspect ratio. This might be a problem for some WMS servers
#    presenting data in a highly inappropriate projection. Do you really
#    want their data?

#%module
#% description: Produces tilings of the source projection for use in the destination region and projection.
#% keywords: raster
#% keywords: tiling
#%end
#%flag
#% key: g
#% description: Produces shell script output
#%end
#%flag
#% key: w
#% description: Produces web map server query string output
#%end
#%option
#% key: region
#% type: string
#% description: Name of region to use instead of current region for bounds and resolution
#%end
#%option
#% key: sourceproj
#% type: string
#% description: Source projection
#% required : yes
#%end
#%option
#% key: sourcescale
#% type: string
#% description: Conversion factor from units to meters in source projection
#% answer : 1
#%end
#%option
#% key: destproj
#% type: string
#% description: Destination projection, defaults to this location's projection
#% required : no
#%end
#%option
#% key: destscale
#% type: string
#% description: Conversion factor from units to meters in source projection
#% required : no
#%end
#%option
#% key: maxcols
#% type: integer
#% description: Maximum number of columns for a tile in the source projection
#% answer: 1024
#%end
#%option
#% key: maxrows
#% type: integer
#% description: Maximum number of rows for a tile in the source projection
#% answer: 1024
#%end
#%option
#% key: overlap
#% type: integer
#% description: Number of cells tiles should overlap in each direction
#% answer: 0
#%end
#%option G_OPT_F_SEP
#% description: Output field separator
#%end

# Data structures used in this program:
# A bounding box:
# 0 -> left, 1-> bottom, 2->right, 3-> top
# A border:
# An array of points indexed by 0 for "x" and 4 for "y" + by number 0, 1, 2, and 3
# A reprojector [0] is name of source projection, [1] is name of destination
# A projection - [0] is proj.4 text, [1] is scale

import sys
import subprocess
import tempfile
import math

from grass.script.utils import separator
from grass.script import core as grass

def bboxToPoints(bbox):
    """Make points that are the corners of a bounding box"""
    points = []
    points.append((bbox['w'], bbox['s']))
    points.append((bbox['w'], bbox['n']))
    points.append((bbox['e'], bbox['n']))
    points.append((bbox['e'], bbox['s']))
    
    return points

def pointsToBbox(points):
    bbox = {}
    min_x = min_y = max_x = max_y = None
    for point in points:
        if not min_x:
            min_x = max_x = point[0]
        if not min_y:
            min_y = max_y = point[1]
        
        if min_x > point[0]:
            min_x = point[0]
        if max_x < point[0]:
            max_x = point[0]
        if min_y > point[1]:
            min_y = point[1]
        if max_y < point[1]:
            max_y = point[1]
    
    bbox['n'] = max_y
    bbox['s'] = min_y
    bbox['w'] = min_x
    bbox['e'] = max_x
    
    return bbox

def project(file, source, dest):
    """Projects point (x, y) using projector"""
    points = []
    ret = grass.read_command('m.proj',
                             quiet = True,
                             flags = 'd',
                             proj_in = source['proj'],
                             proj_out = dest['proj'],
                             sep = ';',
                             input = file)
    
    if not ret:
        grass.fatal(cs2cs + ' failed')
    
    for line in ret.splitlines():
        p_x2, p_y2, p_z2 = map(float, line.split(';'))
        points.append((p_x2 / dest['scale'], p_y2 / dest['scale']))
    
    return points
    
def projectPoints(points, source, dest):
    """Projects a list of points"""
    dest_points = []
    
    input = tempfile.NamedTemporaryFile(mode="wt")
    for point in points:
        input.file.write('%f;%f\n' % \
                             (point[0] * source['scale'],
                              point[1] * source['scale']))
    input.file.flush()
    
    dest_points = project(input.name, source, dest)
    
    return dest_points

def sideLengths(points, xmetric, ymetric):
    """Find the length of sides of a set of points from one to the next"""
    ret = []
    for i in range(len(points)):
        x1, y1 = points[i]
        j = i + 1
        if j >= len(points):
            j = 0
        sl_x = (points[j][0] - points[i][0]) * xmetric
        sl_y = (points[j][1] - points[i][1]) * ymetric
        sl_d = math.sqrt(sl_x * sl_x + sl_y * sl_y)
        ret.append(sl_d)
    
    return { 'x' : (ret[1], ret[3]),
             'y' : (ret[0], ret[2]) }

def bboxesIntersect(bbox_1, bbox_2):
    """Determine if two bounding boxes intersect"""
    bi_a1 = (bbox_1['w'], bbox_1['s'])
    bi_a2 = (bbox_1['e'], bbox_1['n'])
    bi_b1 = (bbox_2['w'], bbox_2['s'])
    bi_b2 = (bbox_2['e'], bbox_2['n'])
    cin = [False, False]
    for i in (0, 1):
        if (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b1[i]) or \
               (bi_a1[i] <= bi_b1[i] and bi_a2[i] >= bi_b2[i]) or \
               (bi_b1[i] <= bi_a1[i] and bi_b2[i] >= bi_a1[i]) or \
               (bi_b1[i] <= bi_a1[i] and bi_b2[i] >= bi_a2[i]):
            cin[i] = True
    
    if cin[0] and cin[1]:
        return True
    
    return False

def main():
    # Take into account those extra pixels we'll be a addin'
    max_cols = int(options['maxcols']) - int(options['overlap'])
    max_rows = int(options['maxrows']) - int(options['overlap'])
    
    # destination projection
    if not options['destproj']:
        dest_proj = grass.read_command('g.proj',
                                       quiet = True,
                                       flags = 'jf').rstrip('\n')
        if not dest_proj:
            grass.fatal(_('g.proj failed'))
    else:
        dest_proj = options['destproj']
    grass.debug("Getting destination projection -> '%s'" % dest_proj)
        
    # projection scale
    if not options['destscale']:
        ret = grass.parse_command('g.proj',
                                 quiet = True,
                                 flags = 'j')
        if not ret:
            grass.fatal(_('g.proj failed'))

        if '+to_meter' in ret:
            dest_scale = ret['+to_meter'].strip()
        else:
            grass.warning(_("Scale (%s) not found, assuming '1'") % '+to_meter')
            dest_scale = '1'
    else:
        dest_scale = options['destscale']
    grass.debug('Getting destination projection scale -> %s' % dest_scale)
    
    # set up the projections
    srs_source = { 'proj' : options['sourceproj'],
                   'scale' : float(options['sourcescale']) }
    srs_dest   = { 'proj' : dest_proj,
                   'scale' : float(dest_scale) }   
    
    if options['region']:
        grass.run_command('g.region',
                          quiet = True,
                          region = options['region'])
    dest_bbox = grass.region()
    grass.debug('Getting destination region')

    # output field separator
    fs = separator(options['separator'])

    # project the destination region into the source:
    grass.verbose('Projecting destination region into source...')
    dest_bbox_points = bboxToPoints(dest_bbox)
    
    dest_bbox_source_points = projectPoints(dest_bbox_points,
                                            source = srs_dest,
                                            dest = srs_source)
    
    source_bbox = pointsToBbox(dest_bbox_source_points)
    
    grass.verbose('Projecting source bounding box into destination...')
    
    source_bbox_points = bboxToPoints(source_bbox)
    
    source_bbox_dest_points = projectPoints(source_bbox_points,
                                            source = srs_source,
                                            dest = srs_dest)
    
    x_metric = 1 / dest_bbox['ewres']
    y_metric = 1 / dest_bbox['nsres']
    
    grass.verbose('Computing length of sides of source bounding box...')
    
    source_bbox_dest_lengths = sideLengths(source_bbox_dest_points,
                                           x_metric, y_metric)
    
    # Find the skewedness of the two directions.
    # Define it to be greater than one
    # In the direction (x or y) in which the world is least skewed (ie north south in lat long)
    # Divide the world into strips. These strips are as big as possible contrained by max_
    # In the other direction do the same thing.
    # Theres some recomputation of the size of the world that's got to come in here somewhere.
    
    # For now, however, we are going to go ahead and request more data than is necessary.
    # For small regions far from the critical areas of projections this makes very little difference
    # in the amount of data gotten.
    # We can make this efficient for big regions or regions near critical points later.
    
    bigger = []
    bigger.append(max(source_bbox_dest_lengths['x']))
    bigger.append(max(source_bbox_dest_lengths['y']))
    maxdim = (max_cols, max_rows)
    
    # Compute the number and size of tiles to use in each direction
    # I'm making fairly even sized tiles
    # They differer from each other in height and width only by one cell
    # I'm going to make the numbers all simpler and add this extra cell to
    # every tile.
    
    grass.message(_('Computing tiling...'))
    tiles = [-1, -1]
    tile_base_size = [-1, -1]
    tiles_extra_1 = [-1, -1]
    tile_size = [-1, -1]
    tileset_size = [-1, -1]
    tile_size_overlap = [-1, -1]
    for i in range(len(bigger)):
	# make these into integers.
	# round up
	bigger[i] = int(bigger[i] + 1)
	tiles[i] = int((bigger[i] / maxdim[i]) + 1)
        tile_size[i] = tile_base_size[i] = int(bigger[i] / tiles[i])
	tiles_extra_1[i] = int(bigger[i] % tiles[i])
	# This is adding the extra pixel (remainder) to all of the tiles:
        if tiles_extra_1[i] > 0:
            tile_size[i] = tile_base_size[i] + 1
        tileset_size[i] = int(tile_size[i] * tiles[i])
	# Add overlap to tiles (doesn't effect tileset_size
        tile_size_overlap[i] = tile_size[i] + int(options['overlap'])
    
    grass.verbose("There will be %d by %d tiles each %d by %d cells" % \
                  (tiles[0], tiles[1], tile_size[0], tile_size[1]))
    
    ximax = tiles[0]
    yimax = tiles[1]
    
    min_x = source_bbox['w']
    min_y = source_bbox['s']
    max_x = source_bbox['e']
    max_y = source_bbox['n']
    span_x = (max_x - min_x)
    span_y = (max_y - min_y)
    
    xi = 0
    tile_bbox = { 'w' : -1, 's': -1, 'e' : -1, 'n' : -1 }
    while xi < ximax:
        tile_bbox['w'] = float(min_x) + (float(xi) * float(tile_size[0]) / float(tileset_size[0])) * float(span_x)
        tile_bbox['e'] = float(min_x) + (float(xi + 1) * float(tile_size_overlap[0]) / float(tileset_size[0])) * float(span_x)
	yi = 0
        while yi < yimax:
            tile_bbox['s'] = float(min_y) + (float(yi) * float(tile_size[1]) / float(tileset_size[1])) * float(span_y)
            tile_bbox['n'] = float(min_y) + (float(yi + 1) * float(tile_size_overlap[1]) / float(tileset_size[1])) * float(span_y)
            tile_bbox_points = bboxToPoints(tile_bbox)
            tile_dest_bbox_points = projectPoints(tile_bbox_points,
                                                  source = srs_source,
                                                  dest = srs_dest)
            tile_dest_bbox = pointsToBbox(tile_dest_bbox_points)
            if bboxesIntersect(tile_dest_bbox, dest_bbox):
                if flags['w']:
                    print "bbox=%s,%s,%s,%s&width=%s&height=%s" % \
                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'], tile_bbox['n'],
                           tile_size_overlap[0], tile_size_overlap[1])
                elif flags['g']:
                    print "w=%s;s=%s;e=%s;n=%s;cols=%s;rows=%s" % \
                          (tile_bbox['w'], tile_bbox['s'], tile_bbox['e'], tile_bbox['n'],
                           tile_size_overlap[0], tile_size_overlap[1])
                else:
                    print "%s%s%s%s%s%s%s%s%s%s%s" % \
                          (tile_bbox['w'], fs, tile_bbox['s'], fs, tile_bbox['e'], fs, tile_bbox['n'], fs,
                           tile_size_overlap[0], fs, tile_size_overlap[1])
            yi += 1
        xi += 1
    
if __name__ == "__main__":
    cs2cs = 'cs2cs'
    options, flags = grass.parser()
    sys.exit(main())