Skip to content

Commit 87ae031

Browse files
author
Rusty Holleman
committed
Rounding bug and openfoam example notebook
1 parent 3b5ea2c commit 87ae031

File tree

2 files changed

+339
-5
lines changed

2 files changed

+339
-5
lines changed

examples/openfoam.ipynb

Lines changed: 324 additions & 0 deletions
Large diffs are not rendered by default.

stompy/spatial/field.py

Lines changed: 15 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -636,7 +636,7 @@ def rectify(self,dx=None,dy=None):
636636
F=newF,projection=self.projection())
637637

638638
def to_grid(self,nx=2000,ny=2000,interp='linear',bounds=None,dx=None,dy=None,
639-
aspect=1.0,max_radius=None,clean=False):
639+
aspect=1.0,max_radius=None,clean=False,strict_dx=False):
640640
""" use the delaunay based griddata() to interpolate this field onto
641641
a rectilinear grid. In theory interp='linear' would give bilinear
642642
interpolation, but it tends to complain about grid spacing, so best to stick
@@ -649,6 +649,8 @@ def to_grid(self,nx=2000,ny=2000,interp='linear',bounds=None,dx=None,dy=None,
649649
interp='qhull': use scipy's delaunay/qhull interface. this can
650650
additionally accept a radius which limits the output to triangles
651651
with a smaller circumradius.
652+
653+
strict_dx: adjust bounds to match dx/dy
652654
"""
653655
if bounds is None:
654656
xmin,xmax,ymin,ymax = self.bounds()
@@ -665,12 +667,20 @@ def to_grid(self,nx=2000,ny=2000,interp='linear',bounds=None,dx=None,dy=None,
665667
# xmin = xmin - (xmin%dx)
666668
# ymin = ymin - (ymin%dy)
667669

670+
# used to truncate and always add 1, but this seems
671+
# better since it does the right thing on an exact
672+
# result.
673+
nx=int( np.ceil((xmax - xmin)/dx) )
674+
ny=int( np.ceil((ymax - ymin)/dy) )
675+
# Previous approach:
668676
# The 1+, -1, stuff feels a bit sketch. But this is how
669677
# CompositeField calculates sizes
670-
nx = 1 + int( (xmax-xmin)/dx )
671-
ny = 1 + int( (ymax-ymin)/dy )
672-
xmax = xmin + (nx-1)*dx
673-
ymax = ymin + (ny-1)*dy
678+
# nx = 1 + int( (xmax-xmin)/dx )
679+
# ny = 1 + int( (ymax-ymin)/dy )
680+
if strict_dx:
681+
xmax = xmin + (nx-1)*dx
682+
ymax = ymin + (ny-1)*dy
683+
# otherwise the user doesn't quite get the requested dx and dy
674684

675685
# hopefully this is more compatible between versions, also exposes more of what's
676686
# going on

0 commit comments

Comments
 (0)