Plot coordinates of a tif image

Hello, I have say 'image.tif' which contains a greyscale information and a .json file containing metadata associated.

I can plot the image using:

[A,R] = readgeoraster('image.tif');
figure
mapshow(A,R)

However, it does come with X or Y axes in a unit I do not understand, value with x10^6, say 600,000 to 605,000 for Y axis. Same fo X. See figure attached from mapshow help page:

https://www.mathworks.com/help/examples/map/win64/OverrideDefaultPropertiesOfTheLineExample_01.png

I would like to transform it into proper longitude/latitude, so °N/°E when I plot the image. In the metadata, there are different projection informationn such as projectedCRS, polygon coordinates and projection epsg/shape/transform values.

I tried using the function worldGrid and geoshow functions, but 1st) it takes ages to interpolate lat/lon and 2nd) my computer does not have enough memory to plot the results using geoshow(lat,lon,A). My images are pretty small but high resolution (20,000 x 20,000 pixels).

So my question, how can I convert my Y and X axes into latitude/longitude degrees on the plot? Any function that would do it? In the end I just need to plot a grayscale image with coordinate.

[Dyuman Joshi]

Try specifying the 'CoordinateSystemType' as 'geographic'.

https://in.mathworks.com/help/map/ref/readgeoraster.html#mw_d0892e98-1c13-4d44-bacc-470c6074ab14

[Madmad]

Error using readgeoraster
The specified coordinate system type, 'geographic', is inconsistent with the projection information associated with
'image.tif'.

Sadly this is not the solution. My data really come with X/Y values not being lat/lon, even if specific projections are available within R.

[Dyuman Joshi]

Could you share the image here? Use the paperclip button to attach.

[Madmad]

I could share but they are quite big (600 MB). I got it from public SAR dataset so you would need to make a free account (https://platform.capellaspace.com/user/login/?redirectConsole=%2Fsearch%2F)

Answer 1

[Cris LaPierre]

Note that there is a projection associated with your SAR data. You need to know this projection in order to convert your scales to lat/lon. This should be incorporated into your geotiff file.

Follow this example to either display your projected data, or to unproject the data.

Project and Display Raster Data

[Madmad]

Yes this is exactly what I have done and said I have done in my initial message

[Z,R] = readgeoraster('image.tif');
p = R.ProjectedCRS;
[x,y] = worldGrid(R);
[lat,lon] = projinv(p,x,y);
figure
geoshow(lat,lon,Z)

The geoshow line fails, because I need somehow more memory, lat, lon, Z are 20,000 x 20,000 approximately.

[Cris LaPierre]

Not exactly the same, as you said you were using mapshow.

Help us out by attaching your image and all the relevant code. You may need to zip the image in order to attach it using the paperclip icon..

[Madmad]

Not exactly the same, as you said you were using mapshow. Sure I initially mentioned mapshow, but I also said that I tried geoshow in my message

• I tried using the function worldGrid and geoshow functions, but 1st) it takes ages to interpolate lat/lon and 2nd) my computer does not have enough memory to plot the results using geoshow(lat,lon,A). My images are pretty small but high resolution (20,000 x 20,000 pixels).

Help us out by attaching your image and all the relevant code. You may need to zip the image in order to attach it using the paperclip icon.

I would love to but my file is too big to be shared. However an easy way to reproduce it is:

lat_t = linspace(18.1626,18.2351,20000);
lon_t = linspace(-73.7939,-73.7183,20000);
[lat,lon] = meshgrid(lat_t,lon_t); % This is similar to my lat/lon after
%  I processed projinv from the worldGrid x and y values
% p = R.ProjectedCRS;
% [x,y] = worldGrid(R);
% [lat,lon] = projinv(p,x,y);
Z = uint16(mapminmax(rand(20000,20000),0,65256)); % Random but same dimension as my real Z
geoshow(lat,lon,Z)
% Out of memory.

I have a memory issue after a few minutes waiting for geoshow (16GB RAM)

In reallity my lat and lon are not consistent along one dimension, so maybe even more complex.

Any comment on how to allow my computer process it would be great!

[Cris LaPierre]

It appears the actual issue is working with large images. Short of increasing your hardware, there is really only one solution - reduce the size of your image.

See here: MATLAB closes when displaying large TIFF image with "geoshow"

[Madmad]

You are right. Two last questions:

1) How much GB do I need to process such figure? Can I calculate it? Because it says "not enough memory", but do not really says how much more I need? Moreover, why can mapshow create a figure in 1 second, but geoshow cannot work on my low spec computer, as I am basically creating the same thing just slightly reprojected?

2) Is there a workaround to get the axes from geoshow (all I need is the lat/lon coordinates of the figure axes), and superpose a mapshow version on top of it (plot NaN with geoshow for example). Or maybe a unit conversion (uint16 for Z and double for lat/lon can be reduced?). Because in the end, what I need is a tif with coherent coordinates, the projection is not really important, I just need an idea of where an area of the image is with a brief look at the x/y axes and grid.

[Cris LaPierre]

1. You whould be able to determine the size of your image using any file browser. In order to load the image into memory, you will need at least that much RAM.
2. You could use the same commands as in the example I shared above, but then plot using mapshow instead of geoshow. The main difference I discovered is that I had to plot lon as X and lat as y.

p = R.ProjectedCRS;
[x,y] = worldGrid(R);
[lat,lon] = projinv(p,x,y);
figure
mapshow(lon,lat,Z)

xlabel('Longitude (degrees)')
ylabel('Latitude (degrees)')

[Cris LaPierre]

I have more details for you.

mapshow is faster because it assumes your image is in cartesian coordinates. The one we have been using - boston.tif - is, so this works. However, it is more likely than not that you geotiff is not. This is why you should use geoshow instead.

geoshow projects the image using the projection info in the geotiff. This can take some time, especially for very large files.