|
|
@@ -108,12 +108,13 @@
|
|
|
# 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
|
|
|
+from grass.exceptions import CalledModuleError
|
|
|
+
|
|
|
|
|
|
def bboxToPoints(bbox):
|
|
|
"""Make points that are the corners of a bounding box"""
|
|
|
@@ -122,9 +123,10 @@ def bboxToPoints(bbox):
|
|
|
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
|
|
|
@@ -133,7 +135,7 @@ def pointsToBbox(points):
|
|
|
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]:
|
|
|
@@ -142,48 +144,57 @@ def pointsToBbox(points):
|
|
|
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"""
|
|
|
+ errors = 0
|
|
|
points = []
|
|
|
- ret = grass.read_command('m.proj',
|
|
|
- quiet = True,
|
|
|
- flags = 'd',
|
|
|
- proj_in = source['proj'],
|
|
|
- proj_out = dest['proj'],
|
|
|
- sep = ';',
|
|
|
- input = file)
|
|
|
-
|
|
|
+ try:
|
|
|
+ ret = grass.read_command('m.proj',
|
|
|
+ quiet=True,
|
|
|
+ flags='d',
|
|
|
+ proj_in=source['proj'],
|
|
|
+ proj_out=dest['proj'],
|
|
|
+ sep=';',
|
|
|
+ input=file)
|
|
|
+ except CalledModuleError:
|
|
|
+ grass.fatal(cs2cs + ' failed')
|
|
|
+
|
|
|
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
|
|
|
-
|
|
|
+ if "*" in line:
|
|
|
+ errors += 1
|
|
|
+ else:
|
|
|
+ p_x2, p_y2, p_z2 = map(float, line.split(';'))
|
|
|
+ points.append((p_x2 / dest['scale'], p_y2 / dest['scale']))
|
|
|
+
|
|
|
+ return points, errors
|
|
|
+
|
|
|
+
|
|
|
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.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
|
|
|
+
|
|
|
+ dest_points, errors = project(input.name, source, dest)
|
|
|
+
|
|
|
+ return dest_points, errors
|
|
|
+
|
|
|
|
|
|
def sideLengths(points, xmetric, ymetric):
|
|
|
"""Find the length of sides of a set of points from one to the next"""
|
|
|
@@ -197,9 +208,9 @@ def sideLengths(points, xmetric, ymetric):
|
|
|
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]) }
|
|
|
+
|
|
|
+ return {'x': (ret[1], ret[3]), 'y': (ret[0], ret[2])}
|
|
|
+
|
|
|
|
|
|
def bboxesIntersect(bbox_1, bbox_2):
|
|
|
"""Determine if two bounding boxes intersect"""
|
|
|
@@ -210,37 +221,46 @@ def bboxesIntersect(bbox_1, bbox_2):
|
|
|
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]):
|
|
|
+ (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'])
|
|
|
-
|
|
|
+
|
|
|
+ if max_cols == 0:
|
|
|
+ grass.fatal(_("It is not possibile to set 'maxcols=%s' and "
|
|
|
+ "'overlap=%s'. Please set maxcols>overlap" %
|
|
|
+ (options['maxcols'], options['overlap'])))
|
|
|
+ elif max_rows == 0:
|
|
|
+ grass.fatal(_("It is not possibile to set 'maxrows=%s' and "
|
|
|
+ "'overlap=%s'. Please set maxrows>overlap" %
|
|
|
+ (options['maxrows'], options['overlap'])))
|
|
|
# destination projection
|
|
|
if not options['destproj']:
|
|
|
dest_proj = grass.read_command('g.proj',
|
|
|
- quiet = True,
|
|
|
- flags = 'jf').rstrip('\n')
|
|
|
+ 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')
|
|
|
+ quiet=True,
|
|
|
+ flags='j')
|
|
|
if not ret:
|
|
|
grass.fatal(_('g.proj failed'))
|
|
|
|
|
|
@@ -252,17 +272,16 @@ def main():
|
|
|
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) }
|
|
|
-
|
|
|
+ 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'])
|
|
|
+ quiet=True,
|
|
|
+ region=options['region'])
|
|
|
dest_bbox = grass.region()
|
|
|
grass.debug('Getting destination region')
|
|
|
|
|
|
@@ -272,52 +291,57 @@ def main():
|
|
|
# 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)
|
|
|
-
|
|
|
+
|
|
|
+ dest_bbox_source_points, errors_dest = projectPoints(dest_bbox_points,
|
|
|
+ source=srs_dest,
|
|
|
+ dest=srs_source)
|
|
|
+
|
|
|
+ if len(dest_bbox_source_points) == 0:
|
|
|
+ grass.fatal(_("There are no tiles available. Probably the output "
|
|
|
+ "projection system it is not compatible with the "
|
|
|
+ "projection of the current location"))
|
|
|
+
|
|
|
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)
|
|
|
-
|
|
|
+
|
|
|
+ source_bbox_dest_points, errors_source = 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]
|
|
|
@@ -326,62 +350,70 @@ def main():
|
|
|
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)
|
|
|
+ # 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:
|
|
|
+ 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
|
|
|
+ # 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" % \
|
|
|
+
|
|
|
+ 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 }
|
|
|
+ tile_bbox = {'w': -1, 's': -1, 'e': -1, 'n': -1}
|
|
|
+
|
|
|
+ if errors_dest > 0:
|
|
|
+ grass.warning(_("During computation %i tiles could not be created" %
|
|
|
+ errors_dest))
|
|
|
+
|
|
|
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
|
|
|
+ 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_points, errors = 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])
|
|
|
+ (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])
|
|
|
+ (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_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()
|
Luca Delucchi