Where do two 2-D arrays begin to overlap each other?

Question

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?

Answer 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)