(Re-)create a matrix of latitude and one of longitude knowing the 4 vertices , taking into account the earth curvature

Dear all,

I wonder if the below can be implemented in MATLAB (or at all).

I have got a satellite image (Z). Z is a matrix of 2000 x 2048. I need to create two matrices, one of longitude and one of latitude, both having the same size of Z (2000 x 2048).

The information I have got to build this matrices are:

Matrix Latitude corners = [47.326; 41.832; 29.8886; 25.460]

Matrix Longitude corners = [3.769; 39.255; 2.226; 30.709]

Reference Ellipsoid Model = WGS84

I have tried several interpolation approaches, but of course, none of them takes into consideration the Earth curvature. Any help would be appreciated!

Thanks!

[John D'Errico]

What about the "matrix" are you trying to build? Just a set of coordinates of lattitude and longitude inside that quadrilateral region, but essentially where the edges are all great circles on the earth? Do you have the mapping toolbox? I would start there.

https://www.mathworks.com/help/map/great-circles-geodesics-and-rhumb-lines.html

[Simone A.]

Dear John, thanks for getting back to me. I do have the mapping toolbox. I didn't understand exactly what you mean by " but essentially where the edges are all great circles on the earth?". The matrix I'm trying to build should provide, per each "pixel", one value of latitude (the latitude matrix) and one of longitude (in the longitude matrix). When these matrices are used to georeference the data matrix (i.e., the actual satellite image) (Z), they should correctly assign the lat/Lon values corresponding to the specific scene within the image. If I take the vertices of the quadrilateral region and I interpolate them to create a 2000 X 2048 matrix, the interpolation would result in an equally spaced grid. However, being this a satellite image, the earth's curvature (at least) should be taken into account when the interpolation is done. I'm not sure if this is "doable", or if other parameters such as satellite altitude are needed. That's why I am seeking for your help! Thanks once again

Answer 1

[Anshuman]

You can create the latitude and longitude matrices using the above information by using MATLAB mapping Toolbox to convert between different coordinate systems and projections.

Here is an example code snippet:

% Load the Mapping Toolbox
% i am assuming you already have the Mapping toolbox or you can install it
% by typing "matlab.addons.install" in the cmd window of MATLAB

% Define the image size
imageSize = [2000, 2048];

% Define the latitude and longitude corners in geographic coordinates
latitudeCorners = [47.326; 41.832; 29.8886; 25.460];
longitudeCorners = [3.769; 39.255; 2.226; 30.709];

% Define the reference ellipsoid model
refEllipsoid = wgs84Ellipsoid();

% Define the projection for the satellite image
satelliteProjection = mapProjection('Equirectangular', 'geoid', refEllipsoid);

% Define the X and Y limits for the image
xLimits = [1, imageSize(2)];
yLimits = [1, imageSize(1)];

% Convert the latitude and longitude corners to projected coordinates
[xCorners, yCorners] = project(satelliteProjection, latitudeCorners, longitudeCorners);

% Create grids of X and Y coordinates
[X, Y] = meshgrid(linspace(xLimits(1), xLimits(2), imageSize(2)), ...
                  linspace(yLimits(1), yLimits(2), imageSize(1)));

% Convert the X and Y coordinates to latitude and longitude
[lat, lon] = inverse(satelliteProjection, X, Y);

% Display the latitude and longitude matrices
disp(lat);
disp(lon);

The project and inverse functions are used to convert between geographic coordinates and projected coordinates, taking into account the curvature of the earth.

Hope it helps!

[Simone A.]

Hi @Anshuman, thanks for getting back to me. Although I have got the mapping Toolbox installed, I get the following error. Can you help with it? I cannot find the "mapProjection" function anywhere.

>> satelliteProjection = mapProjection('Equirectangular', 'geoid', refEllipsoid);

Unrecognized function or variable 'mapProjection'.

[Anshuman]

You can also try to use the projectedCRS function in place of mapProjection:

% Create a projected coordinate reference system object for an equirectangular projection:
satelliteProjection = projectedCRS('Equirectangular', 'geographicCRS', geocrs('WGS 84'), 'geoid', refEllipsoid);

% Create a geographic coordinate reference object for the corners of the image:
imageGeographicCRS = geocrs('WGS 84', 'Degrees');

% Project the corners of the image to the equirectangular projection:
[xCorners, yCorners] = project(imageGeographicCRS, satelliteProjection, longitudeCorners, latitudeCorners);

[Simone A.]

Hi, thanks for getting back. It is still not working:

I first get:

% Create a projected coordinate reference system object for an equirectangular projection:
>> satelliteProjection = projectedCRS('Equirectangular', 'geographicCRS', geocrs('WGS 84'), 'geoid', refEllipsoid);
Error using geocrs
Argument must specify a geographic coordinate reference system.

So I changed it by adding:

>> satelliteProjection = projectedCRS('Equirectangular', 'geographicCRS', geocrs(4326,'Authority','EPSG'), 'geoid', refEllipsoid);

But I get the error:

% Unrecognized function or variable 'projectedCRS'.