Join geo-referenced matrices in MATLAB

Question

I have two geo-referenced matrices (gridded climate datasets) with different spatial extent and resolution (i.e., spatial coverage of one pixel on the surface of the Earth) that I would like to join. Let's call them referenceMatrix and targetMatrix. The matrices can be loaded into MATLAB as geotiffs or just matrices with corresponding grids of latitudes/longitudes per pixel.

What I want is to fill NaN-pixels in the targetMatrix with the corresponding value from the referenceMatrix. I can do it using a for-loop that looks through the pixels one by one, and fills the NaNs in the targetMatrix with data from the referenceMatrix based on the nearest pixel. The method I am using to localize the nearest pixel in space is described here:
How to find the nearest points to given coordinates with MATLAB?

However I need to do this with thousands of matrices, and the for-loop is too slow. Therefore I would like to use logical indexing, e.g.

targetMatrix(isnan(targetMatrix)) = referenceMatrix(isnan(targetMatrix))

with the added ability to match pixels in the matrices based on their latitude and longitude. Can any of you point me in a direction with an example of comparing matrices with different extent based on their georeference?

An example of input data and the wanted output below

targetMatrix = [1,  NaN, 3; 
                NaN, 5,  6]; 
referenceMatrix = [10, 20, 30, 40; 
                   50, 60, 70, 80]; 
referenceLatitude = [13.3, 13.3, 13.3, 13.3; 
                     14.1, 14.1, 14.1, 14.1]; 
referenceLongitude = [3.2, 4.2, 5.2, 6.2; 
                      3.2, 4.2, 5.2, 6.2];  
targetLatitude = [13.4, 13.4, 13.4; 
                  13.9, 13.9, 13.9]; 
targetLongitude = [3.1, 3.6, 4.1; 
                   3.1, 3.6, 4.1]; 

wantedOutput = [ 1, 10, 3; 
                50,  5, 6];

The wanted output consists of the original values from the targetMatrix where the NaNs are filled with the nearest (in space) value from the referenceMatrix, i.e., 10 and 50.

The NaN in the targetMatrix(1,2) has the lat/lon coordinates (13.4, 3.6). The value 10 in the referenceMatrix(1,1) has the lat/lon coordinates (13.3, 3.2) and is thus the point located closest in space to targetMatrix(1,2). The value 20 in the referenceMatrix has the lat/lon coordinates (13.3, 4.2) and the point is therefore further away in space from the targetMatrix(1,2). — user3275421

Sardar Usama Sardar Usama · Accepted Answer · 2018-06-12T10:26:32

Find which entries of targetMatrix are to be replaced using isnan. Convert the lat/lon from degrees to radians and then to cartesian coordinates using sph2cart to get actual geodesic distances. Then use knnsearch^⬥ to find the indices of the points from the reference coordinates which are nearest to relevant targeted coordinates. Use these indices to extract the relevant entries from referenceMatrix and replace NaNs with them.

nanent = isnan(targetMatrix);
[tarX, tarY] = sph2cart(targetLongitude*pi/180, targetLatitude*pi/180, 1);
[refX, refY] = sph2cart(referenceLongitude*pi/180, referenceLatitude*pi/180, 1);
tmptar = [tarY(nanent) tarX(nanent)];
tmpref = [refY(:) refX(:)];
ind = knnsearch(tmpref, tmptar);
wantedOutput = targetMatrix;
wantedOutput(nanent) = referenceMatrix(ind);

In this case, converting lat/lon to cartesian coordinates prior to the use of knnsearch is tested^ⵜ to sped up knnsearch than otherwise.

knnsearch finds euclidean distance by default. You can also do this using a combination of pdist2^⬥ and min instead of knnsearch (in line 6). i.e.

[~, ind] = min(pdist2(tmptar, tmpref), [], 2);

or you can use desearchn in line 6. i.e.

ind = dsearchn(tmpref, tmptar);

But knnsearch is tested^ⵜ to be faster than dsearchn in this case.

_{^⬥knnsearch and pdist2 require ≥ R2010a with Stats and ML Toolbox. Use dsearchn if you do not have that.}
_{^ⵜ All tests are done by the OP.}

Join geo-referenced matrices in MATLAB

1 Answers