Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

allow trim of internal_holes (houses) when producing geotiff outputs using Make_Geotif #167

Merged
merged 2 commits into from
Mar 10, 2018
Merged

allow trim of internal_holes (houses) when producing geotiff outputs using Make_Geotif #167

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
7fdf0f0
allow trimming of internal_holes (houses) when producing geotiff outp…
davekennewell Feb 7, 2018
12b8fe8
added testcase where removal of leading spaces is allowed (occurs on …
davekennewell Mar 9, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
8 changes: 5 additions & 3 deletions anuga/geospatial_data/tests/test_geospatial_data.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,10 +33,12 @@ def test_0(self):

# Check __repr__
# FIXME (Ole): Is this really machine independent?
rep = `G`
rep = G.__repr__()
ref = '[[ 1. 2.1]\n [ 3. 5.3]]'
msg = 'Representation %s is not equal to %s' % (rep, ref)
assert rep == ref, msg
# it seems the leading spaces are removed on some conda builds, so lets allow that too:
ref_conda = '[[1. 2.1]\n [3. 5.3]]'
msg = 'Representation %s is not equal to %s or %s' % (rep, ref, ref_conda)
assert ((rep == ref) or (rep == ref_conda)), msg

#Check getter
assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
Expand Down
17 changes: 16 additions & 1 deletion anuga/utilities/plot_utils.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -957,6 +957,7 @@ def Make_Geotif(swwFile=None,
min_allowed_height=1.0e-05,
output_dir='TIFS',
bounding_polygon=None,
internal_holes=None,
verbose=False,
k_nearest_neighbours=3,
creation_options=[]):
Expand All @@ -983,6 +984,7 @@ def Make_Geotif(swwFile=None,
min_allowed_height -- Minimum allowed height from ANUGA
output_dir -- Write outputs to this directory
bounding_polygon -- polygon (e.g. from read_polygon) If present, only set values of raster cells inside the bounding_polygon
internal_holes -- a list of polygons. If present, do not set values of raster cells inside these polygons.
k_nearest_neighbours -- how many neighbours to use in interpolation. If k>1, inverse-distance-weighted interpolation is used
creation_options -- list of tif creation options for gdal, e.g. ['COMPRESS=DEFLATE']
"""
Expand Down Expand Up @@ -1104,11 +1106,19 @@ def myInterpFun(quantity):
num += quantity[nn_inds[:,i]]*nn_wts[:,i]
return (num/denom)

if(bounding_polygon is not None):
if bounding_polygon is not None:
# Find points to exclude (i.e. outside the bounding polygon)
from anuga.geometry.polygon import outside_polygon
cut_points = outside_polygon(gridXY_array, bounding_polygon)

hole_points_list = []
if internal_holes is not None:
# Find points to exclude (i.e. inside the internal_holes)
from anuga.geometry.polygon import inside_polygon
for hole in internal_holes:
cut_holes = inside_polygon(gridXY_array, hole)
hole_points_list.append(cut_holes)

# Loop over all output quantities and produce the output
for myTSindex, myTSi in enumerate(myTimeStep):
if(verbose):
Expand Down Expand Up @@ -1150,6 +1160,11 @@ def myInterpFun(quantity):
# Cut the points outside the bounding polygon
gridq[cut_points] = numpy.nan

if (internal_holes is not None) and (len(hole_points_list[0]) > 0):
# Cut the points inside the hole polygons
for hole_points in hole_points_list:
gridq[hole_points] = numpy.nan

# Make name for output file
if(myTS != 'pointData'):
output_name = output_dir + '/' +\
Expand Down