Skip to content

added masking of erroneous values #13

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

Merged
pp-mo merged 3 commits into from
Feb 12, 2015
Merged

added masking of erroneous values #13

Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
f9359f3
added masking of erroneous values
esc24 Feb 6, 2015
706cacc
review actions
esc24 Feb 11, 2015
6588dd5
increased tolerance and added testing
esc24 Feb 11, 2015
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 56 additions & 0 deletions lib/iris/analysis/cartography.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -858,6 +858,43 @@ def _transform_distance_vectors(src_crs, x, y, u_dist, v_dist, tgt_crs):
return u2_dist, v2_dist


def _transform_distance_vectors_tolerance_mask(src_crs, x, y, tgt_crs):
"""
Return a mask that can be applied to data array to mask elements
where the magnitude of vectors are not preserved due to numerical
errors introduced by the tranformation between coordinate systems.

Args:
* src_crs (`cartopy.crs.Projection`):
The source coordinate reference systems.
* x, y (array):
Locations of each vector defined in 'src_crs'.
* tgt_crs (`cartopy.crs.Projection`):
The target coordinate reference systems.

Returns:
2d boolean array that is the same shape as x and y.

"""
if x.shape != y.shape:
raise ValueError('Arrays do not have matching shapes. '
'x.shape is {}, y.shape is {}.'.format(x.shape,
y.shape))
ones = np.ones(x.shape)
zeros = np.zeros(x.shape)
u_one_t, v_zero_t = _transform_distance_vectors(src_crs, x, y,
ones, zeros, tgt_crs)
u_zero_t, v_one_t = _transform_distance_vectors(src_crs, x, y,
zeros, ones, tgt_crs)
# Squared magnitudes should be equal to one within acceptable tolerance.
sqmag_1_0 = u_one_t**2 + v_zero_t**2

@pp-mo pp-mo Feb 10, 2015

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Only a tiny point, but I think this would all look simpler with some renames, so that these 2 lines would read (something like):

sqmag_10 = u_10_t**2 + v_10_t**2
sqmag_01 = u_01_t**2 + v_01_t**2

sqmag_0_1 = u_zero_t**2 + v_one_t**2
mask = np.logical_not(

@pp-mo pp-mo Feb 10, 2015

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think at this point we may also need to allow for possible NaNs in the transformed vectors.

@esc24 esc24 Feb 10, 2015

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fortunately np.isclose takes care of these for you.

np.logical_and(np.isclose(sqmag_1_0, ones, atol=1e-3),
np.isclose(sqmag_0_1, ones, atol=1e-3)))
return mask


def rotate_winds(u_cube, v_cube, target_cs):
"""
Transform wind vectors to a different coordinate system.
Expand Down Expand Up @@ -975,12 +1012,26 @@ def rotate_winds(u_cube, v_cube, target_cs):
else:
dims = x_dims

# Transpose x, y 2d arrays to match order in cube's data.
if dims[0] > dims[1]:
x = x.transpose()

@pp-mo pp-mo Feb 10, 2015

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not convinced this is correct..
If we transpose x at this point, then below we will get similarly transposed xt and yt.
( cf. https://github.com/esc24/iris/blob/f9359f3593bdaac1c6753ceca65e40efadc43ff9/lib/iris/analysis/cartography.py#L1057

    # Calculate new coords of locations in target coordinate system.
    xyz_tran = target_crs.transform_points(src_crs, x, y)
    xt = xyz_tran[..., 0].reshape(x.shape)
       ...
    xt_coord = iris.coords.AuxCoord(xt,
                                    standard_name='projection_x_coordinate',
                                    coord_system=target_cs)

But when when assign these back into the cube, we just use 'dims' again, which we have not swapped around ...
( cf. https://github.com/esc24/iris/blob/f9359f3593bdaac1c6753ceca65e40efadc43ff9/lib/iris/analysis/cartography.py#L1074

 ut_cube.add_aux_coord(xt_coord, dims)

)

So it seems to me we may not need to make this change at all.
If we do, I guess we need a test to confirm it works properly.

@esc24 esc24 Feb 10, 2015

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good catch @pp-mo. This was a last minute thing I noticed relating to the data, but coords will need looking at too. I'll give it some thought.

@esc24 esc24 Feb 11, 2015

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've given it thought and I'm now happy. The transpose is necessary as mesh grid will always give yx ordering.

y = y.transpose()

# Create resulting cubes.
ut_cube = u_cube.copy()
vt_cube = v_cube.copy()
ut_cube.rename('transformed_{}'.format(u_cube.name()))
vt_cube.rename('transformed_{}'.format(v_cube.name()))

# Calculate mask based on preservation of magnitude
mask = _transform_distance_vectors_tolerance_mask(src_crs, x, y,
target_crs)
apply_mask = mask.any()
if apply_mask:
# Make masked arrays to accept masking.

@pp-mo pp-mo Feb 10, 2015

Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think you can leave this to the end, and apply the mask to the entire result with broadcasting.
Likewise, you don't need to convert the results to masked until the very end, as the transformed vector values that make the individual slices do not have a mask.

However, the transformed results can have NaNs in, and we should also convert those to masked.
N.B. those cases also don't vary between the slices -- i.e. the NaNs will depend only on the X and Y locations, and not the vector data.

@esc24 esc24 Feb 10, 2015

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You coudl do it at the end with broadcasting, but it's tricky to do it in an additive way i.e. not unset any existing mask. I'd love it if you could you show me how to do this.

@esc24 esc24 Feb 11, 2015

Copy link
Copy Markdown
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I looked at this and its messy as simple broadcasting provided by numpy doesn't cut it. You need to take into account the dimension mapping.

ut_cube.data = ma.asanyarray(ut_cube.data)
vt_cube.data = ma.asanyarray(vt_cube.data)

# Project vectors with u, v components one horiz slice at a time and
# insert into the resulting cubes.
shape = list(u_cube.shape)
Expand All @@ -995,6 +1046,11 @@ def rotate_winds(u_cube, v_cube, target_cs):
u = u_cube.data[index]
v = v_cube.data[index]
ut, vt = _transform_distance_vectors(src_crs, x, y, u, v, target_crs)
if apply_mask:
ut = ma.asanyarray(ut)
ut[mask] = ma.masked
vt = ma.asanyarray(vt)
vt[mask] = ma.masked
ut_cube.data[index] = ut
vt_cube.data[index] = vt

Expand Down