Unable to plot raster using meshm

Hi,

I have the following demo code (using load topo) that works perfectly:

InitialAngle = -25;
Number_of_pannels = 3;

figure
for p = 1:Number_of_pannels
  subplot(1,Number_of_pannels,p)
  %-plot coastlines, and select vantage point-%
  load coastlines
  axesm('ortho','origin',[0 ((360/Number_of_pannels)*p)-InitialAngle]); 
  %-pot info-%
  load topo
  hs = meshm(topo,topolegend,size(topo));
  %-make it pretty-%
  axis off;
  gridm on; 
  framem on;
  plotm(coastlat,coastlon)
end

It produces the following figure:

I want to do the same, but instead of using the raster topo (which basically comes from matlab), I want to show my actual data (that comes from NASA in a tif file).

I read the data using the following:

projection_em = geotiffinfo('EM2001.tif');
[em_raster_tiff,em_reference_tiff] = geotiffread('EM2001.tif');

And then I try to plot it:

hs = meshm(em_raster_tiff,em_reference_tiff);

It gives the following error:

Error using internal.map.convertToGeoRasterRef (line 39)

Function MESHGRAT expected input argument 2, 'R', to be a 3-element

referencing vector, a 3-by-2 referencing matrix, or a (scalar)

geographic raster reference object.

Error in meshgrat (line 139)

R = internal.map.convertToGeoRasterRef( ...

Error in meshm (line 111)

[lat,lon] = meshgrat(Z, R, gratsize);

However,I'm using the "R" that I obtain from geotiffread... and it looks like this:

em_reference_tiff =

MapCellsReference with properties:

XWorldLimits: [-19640076.5206259 19730423.4793741]

YWorldLimits: [-9996423.89867413 9298576.10132587]

RasterSize: [38590 78741]

RasterInterpretation: 'cells'

ColumnsStartFrom: 'north'

RowsStartFrom: 'west'

CellExtentInWorldX: 500

CellExtentInWorldY: 500

RasterExtentInWorldX: 39370500

RasterExtentInWorldY: 19295000

XIntrinsicLimits: [0.5 78741.5]

YIntrinsicLimits: [0.5 38590.5]

TransformationType: 'rectilinear'

CoordinateSystemType: 'planar'

I'm clueless... I see that in the "topo" raster that Matlab uses as an example topolegend = [1 90 0], but I don't know what those numbers refer to.

I apreciate any help.

Cheers!

LG

PS: forgot to mention: mapshow(em_raster_tiff,jet,em_reference_tiff) works, but using mapshow I only get a square map with a stardard projection and no the 3D looking views of the Earth that I'm looking for.

Answer 1

[Luis J Gilarranz]

So, not that is mention anywhere in the Matlab documentation, but I found the RefMatrix the documentation is referring to inside the structure projection_em. Therefore, using the following line of code I no longer get that error:

hs = meshm(em_raster_tiff,projection_em.RefMatrix);

In return, I get the following error:

Error using constructGeoRasterReference (line 45)

In combination with the number of rows in the raster grid, 38590, the

input referencing vector or matrix implies latitude limits that extend

outside the interval [-90 90] degrees.

projection_em.RefMatrix =

1.0e+07 *

0 -0.0001

0.0001 0

-1.9640 0.9299

So, yes, the RefMatrix exceeds the -90 90 boundary, but this is what I got after

projection_em = geotiffinfo('EM2001.tif');

Any ideas?

[Luis J Gilarranz]

the structure projection_em seems to have some useful information, but I don't know how to capitalize on it:

projection_em =

struct with fields:

Filename: '/EM2001.tif'

FileModDate: '21-Jan-2020 17:34:36'

FileSize: 248110757

Format: 'tif'

FormatVersion: []

Height: 38590

Width: 78741

BitDepth: 32

ColorType: 'grayscale'

ModelType: 'ModelTypeProjected'

PCS: ''

Projection: ''

MapSys: ''

Zone: []

CTProjection: 'CT_Sinusoidal'

ProjParm: [7×1 double]

ProjParmId: {7×1 cell}

GCS: 'Unknown datum based upon the Authalic Sphere'

Datum: 'Not specified (based on Authalic Sphere)'

Ellipsoid: 'Sphere'

SemiMajor: 6371000

SemiMinor: 6371000

PM: 'Greenwich'

PMLongToGreenwich: 0

UOMLength: 'metre'

UOMLengthInMeters: 1

UOMAngle: 'degree'

UOMAngleInDegrees: 1

TiePoints: [1×1 struct]

PixelScale: [3×1 double]

SpatialRef: [1×1 map.rasterref.MapCellsReference]

RefMatrix: [3×2 double]

BoundingBox: [2×2 double]

CornerCoords: [1×1 struct]

GeoTIFFCodes: [1×1 struct]

GeoTIFFTags: [1×1 struct]

If we look at projection_em.CornerCoords:

X: [-1.9640e+07 1.9730e+07 1.9730e+07 -1.9640e+07]

Y: [9.2986e+06 9.2986e+06 -9.9964e+06 -9.9964e+06]

Row: [0.5000 0.5000 3.8590e+04 3.8590e+04]

Col: [0.5000 7.8742e+04 7.8742e+04 0.5000]

Lat: [83.6241 83.6241 -89.9000 -89.9000]

Lon: [-150.5097 157.8262 145.5996 -40.0659]

so it seems that there's a way to go from the X, Y to Lat Lon