Geotiffread & getting latlon info

Hello,

I am using EO-1 Hyperion data (hyperspectral) as a geotiff image data of Agatti Islands, Lakshadweep, India. I used geotiffread to read the image data.. using [img, RfrncMtrx, BndgBx]=geotiffread(phileName{:}); When I want to convert from the row/col to lat/lon by using pix2latlon, I get the numbers in map-scale. (Actually the values in the reference matrix are also huge!!!???). The BndgBx=[179100, 1155270; 206430, 1244100]..

Kindly let me know where I am going wrong?

Cheers

Raghu

[Luca Brocca]

info = geotiffinfo(namefile);
[x,y]=pixcenters(info);

[Hyunglok Kim]

Luka's answer is the best.

[Nirajan Luintel]

I used to calculate from boundary box. Its a lot easier. Thanks a lot

Answer 1

[Kelly Luetkemeyer]

Here is updated code to adress this request using recent (R2022a) features in Mapping Toolbox:

%% Read image file using readgeoraster
fname = "boston.tif";
[A,R] = readgeoraster(fname);

%% Create grid of X,Y values
[x,y] = worldGrid(R);

%% Convert grid of X,Y values to latitude/longitude
[lat,lon] = projinv(R.ProjectedCRS,x,y);

%% Plot latitude/longitude min/max on geographic axes to confirm
[latlim,lonlim] = geoquadline(lat,lon);
figure
geoplot(latlim([1 2 2 1 1]),lonlim([1,1,2,2,1]),"r","LineWidth",2)
geobasemap satellite
title("geographic axes")
figure
mapshow(A,R)
title("boston.tif")

********************************

Hi Raghu,

As background information, the coordinates of your image are in a projected coordinate system rather than a geographic coordinate system. You can use the function GEOTIFFINFO to return the information structure.

You will see the 'ModelType' field set to 'ModelTypeProjected'; therefore, the RefMatrix is also referenced to map coordinates.

Beginning in R2011a, you will also see a new field, SpatialRef, which in this case will contain a spatialref.MapRasterReference object.

If you want to find the geographic coordinates from the map coordinates, then you need to use the function PROJINV to unproject the coordinates to latitude and longitude.

Here is how you would find the latitude and longitude values for the center of the first pixel in boston.tif:

info = geotiffinfo('boston.tif');
[x,y] = pix2map(info.RefMatrix, 1, 1);
[lat,lon] = projinv(info, x,y)

Beginning in R2011a, you can use:

[x,y] = pix2map(info.SpatialRef, 1, 1)

or

R = info.SpatialRef;
[x,y] = R.intrinsicToWorld(1,1);

followed by:

    [lat,lon] = projinv(info, x,y);

I hope this help you.

-Kelly

[Raghavendra Mupparthy]

Dear Kelly,

Thank you for the reply and the answer. I finally landed up writing an entire function to convert UTM map scale to lat/lon, till I saw that there is projinv. Thank you I learnt a shorter technique..

I have R2010a, but I do not see info.SpatialRef when I get the geotiffinfo of an image. I am having Mapping toolbox version 3.1.. I tried the same with boston.tif.. still the attribute SpatialRef is not to be seen.. I hope I am not missing any updates..

Cheers

Raghu

[Kelly Luetkemeyer]

Hi Raghu,

You need to have R2011a (Mapping Toolbox version 3.3) in order to see the change to GEOTIFFINFO which adds the SpatailRef field, but for your particular case, the RefMatrix is sufficient.

Best wishes,

-Kelly

[Ehsan Jalilvand]

Dear Kelly

Hi

I want to extract lat and long as a vector or array for all pixels (not just a point!). I want to do this in order to plot the geotiff image by surfm function. Actually I try to use other tools like: mapshow(data,R)or geoshow(data,info.RefMatrix)but they don't work either. So I try to add two for loops to your code as below to solve the problem.

clear

clc

[data, R] = geotiffread('MCD43B3.A2015001.h22v05.005.2015027220419_MOD_Grid_BRDF.tif');

info = geotiffinfo('MCD43B3.A2015001.h22v05.005.2015027220419_MOD_Grid_BRDF.tif');

[AY, AX]= size (data);

for i= 1:AY

      for j=1:AX

[y(i),x(j)] = pix2map(info.RefMatrix, i, j);

[lat(i,j),lon(i,j)] = projinv(info, y(i),x(j));

      end

end

but since the matrix size is 1114x754 the code takes forever to run, how can I solve this problem?

P.S: my image is a MODIS albedo product (MCD43B3) that I converted it to a geotiff using HEG toolbox (HDF-EOS to GeoTIFF Conversion Tool).

Thanks

Ehsan

[TT1981]

I know this is a really old post, but since I was struggling with the same as OP and there are not many posts online related to this, I would like to include a solution to shorten computing time:

 [AY, AX]= size (data);
 lon = zeros(size(data));
 lat = zeros(size(data));
 x = zeros(1, AX);
 y = zeros(1, AY);

 height = info.Height; % Integer indicating the height of the image in pixels
 width = info.Width; % Integer indicating the width of the image in pixels
 [rows,cols] = meshgrid(1:height,1:width);
 % Getting the latitude and longitude of the points
 [x,y] = pix2map(info.RefMatrix, rows, cols);
 [lat,lon] = projinv(info, y,x);

This way you do not need to loop through all the points.

[Tom Bell]

You made my day! Thank you!

Answer 2

[Reema Alhassan]

hello, when I'm using projinv(info, y,x); function I'm getting an error says the GeoTIFF structure PROJ can't be used with the functions PROJFWD or PROJINV if you could help me please ..

thank you

[Emily T. Griffiths]

I got to the bottom of this post, only to have the same question as Reema...

[Sophia Barth]

Had anyone found a solution for this problem yet? Thanks in advance.

[Walter Roberson]

Which release are you using?

What happens if you test with

info = geotiffinfo('boston.tif');
[x,y] = pix2map(info.RefMatrix, 1, 1);
[lat,lon] = projinv(info, x,y)

[Sophia Barth]

Thanks for the prompt reply, I use R2020a.

I exported my image from Google Earth Engine and its in a projected coordinate system ('ModelTypeProjected').

When I use this code for my image the following comes up:

Error using proj2gtif (line 17)

The GeoTIFF structure PROJ cannot be used with functions PROJFWD or PROJINV.

Error in projaccess (line 40)

gtif = proj2gtif(proj);

Error in projinv (line 73)

[lat, lon] = projaccess('inv', proj, x, y);

Error in mathworksanswer (line 6)

[lat,lon] = projinv(info, x,y)

[Walter Roberson]

Can you post code to fetch the data so we can be sure we are working with the same structure?

[Sophia Barth]

This is the code i used to export the image:

https://code.earthengine.google.com/c00424acb0f9d46599cf39fd4353da4a

If you have an account for Google Earth Engine it should work when you click on the link. If the link somehow doesnt work this is the code:

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}
//Filter images
var dataset = ee.ImageCollection('COPERNICUS/S2')
.filterBounds(geometry5)
//.filterDate('2016-07-23','2016-07-27')
.filterDate('2020-07-01','2020-07-30')
.filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than',20)
.map(maskS2clouds)
.median()
.select(['B12','B8','B4'])

//visualization parameter
var SWIRVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B12', 'B8', 'B4'],
};

//output image, add to map
var final_image = dataset.clip(geometry5)
Map.addLayer(final_image,SWIRVis,'Image2016-06');
//Map.addLayer(final_image,VisParam,'image201607_1_fire');

 //Export the image, specifying scale and region
Export.image.toDrive({
  image: final_image,
  description: 'image202007_swir_small',
  scale: 20,
  crs: 'EPSG:3857',
  region: geometry5
});

I hope that helps, thanks again!

[Sophia Barth]

As a next step I dowloaded the image from Googe Drive and used the follwoing code as you suggested:

info = geotiffinfo('image202007_swir_small.tif');
[x,y] = pix2map(info.RefMatrix, 1, 1);
[lat,lon] = projinv(info, x,y)

Answer 3

[Ran Zask]

I think it should be [lat,lon] = projinv(info, x,y);

[Qishun Ran]

I think you are right, it should be [lat,lon] = projinv(info, x,y); and

it should be [rows,cols] = meshgrid(1:width,1:height);

Answer 4

[Andres Rey Sanchez]

Some corrections are needed to the answers above. To get the correct matrices of coordinates (lat, lon) from your geotiff the code should be:

info = geotiffinfo('boston.tif');
height = info.Height; % Integer indicating the height of the image in pixels
width = info.Width; % Integer indicating the width of the image in pixels
[cols,rows] = meshgrid(1:width,1:height);
[x,y] = pix2map(info.RefMatrix, rows, cols);
[lat,lon] = projinv(info, x,y);