6

I'm working with model output at the moment, and I can't seem to come up with a nice way of combining two arrays of data. Arrays A and B store different data, and the entries in each correspond to some spatial (x,y) point -- A holds some parameter, and B holds model output. The problem is that B is a spatial subsection of A -- that is, if the model were for the entire world, A would store the parameter at each point on the earth, and B would store the model output only for those points in Africa.

So I need to find how much B is offset from A -- put another way, I need to find the indexes at which they start to overlap. So if A.shape=(1000,1500), is B the (750:850, 200:300) part of that, or the (783:835, 427:440) subsection? I have arrays associated with both A and B which store the (x,y) positions of the gridpoints for each.

This would seem to be a simple problem -- find where the two arrays overlap. And I can solve it with scipy.spatial's KDTree simply enough, but it's very slow. Anyone have any better ideas?

BernardShaw
  • 63
  • 3
  • 1
    The answer to this question really depends on the type of grid you have. Is it an equally spaced Cartesian grid? – Sven Marnach Nov 11 '10 at 09:07
  • 1
    Are the coordinates ascending both horizontally and vertically for the positions of the gridpoints of A and B? If so, then just do some binary searches. – Justin Peel Nov 13 '10 at 03:57
  • If this is a global climate map, does B correspond to a geographic region like Africa (i.e. something that would not be exactly rectangular in shape) or is it just a rectangular subsection? – Jordan Reiter Nov 22 '10 at 19:40

3 Answers3

1

I have arrays associated with both A and B which store the (x,y) positions of the gridpoints for each.

In that case, the answer should be fairly simple...

Are the two grids strictly on the same gridding scheme? Assuming they are, you can just do something like:

np.argwhere((Ax == Bx.min()) & (Ay == By.min()))

Assuming the world coordinates of the two grids increase in the same direction as the indicies of the grids, this gives the lower left corner of the subsetted grid. (And if they don't increase in the same direction (i.e. negative dx or dy), it just gives one of the other corners)

In the example below, we could obviously just calculate the proper indicies from ix = (Bxmin - Axmin) / dx, etc, but assuming you have a more complex gridding system, this will still work. However, this is assuming that the two grids are on the same gridding scheme! It's slightly more complex if they're not...

import numpy as np

# Generate grids of coordinates from a min, max, and spacing
dx, dy = 0.5, 0.5

# For the larger grid...
Axmin, Axmax = -180, 180
Aymin, Aymax = -90, 90

# For the smaller grid...
Bxmin, Bxmax = -5, 10
Bymin, Bymax = 30, 40

# Generate the indicies on a 2D grid
Ax = np.arange(Axmin, Axmax+dx, dx)
Ay = np.arange(Aymin, Aymax+dy, dy)
Ax, Ay = np.meshgrid(Ax, Ay)

Bx = np.arange(Bxmin, Bxmax+dx, dx)
By = np.arange(Bymin, Bymax+dy, dy)
Bx, By = np.meshgrid(Bx, By)

# Find the corner of where the two grids overlap...
ix, iy = np.argwhere((Ax == Bxmin) & (Ay == Bymin))[0]

# Assert that the coordinates are identical.
assert np.all(Ax[ix:ix+Bx.shape[0], iy:iy+Bx.shape[1]] == Bx) 
assert np.all(Ay[ix:ix+Bx.shape[0], iy:iy+Bx.shape[1]] == By)
Joe Kington
  • 275,208
  • 71
  • 604
  • 463
0

Can you say more? What model are you using? What are you modelling? How is it computed?

Can you make the dimensions match to avoid the fit? (i.e. if B doesn't depend on all of A, only plug in the part of A that B models, or compute boring values for the parts of B that wouldn't overlap A and drop those values later)

user500198
  • 111
  • 1
  • 5
  • It's a global climate model. "A" stores how much velocities at each grid cell have been rotated (for whatever reason), and "B" are the actual velocities. The size and locations of the arrays are the same for all my analysis, so I really only need to figure out the overlap once. As I alluded to, I solved it in a sort of brute-force fashion, which works for what I need to do. So at this point it's more of an abstract intellectual exercise: how does one figure out where one grid, being a slice of another, overlaps the original? – BernardShaw Nov 11 '10 at 04:46
0

I need to find the indexes at which they start to overlap

So are you looking for indexes from A or from B? And is B strictly rectangular?

Finding the bounding box or convex hull of B is really cheap.

Neal Ehardt
  • 10,334
  • 9
  • 41
  • 51