Unable to plot raster using meshm

3 次查看（过去 30 天）
Luis J Gilarranz
Luis J Gilarranz 2020-1-21
评论： Luis J Gilarranz 2020-1-21
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.

请先登录，再回答此问题。

回答（1 个）

Luis J Gilarranz
Luis J Gilarranz 2020-1-21
编辑：Luis J Gilarranz 2020-1-21
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?
  1 个评论
Luis J Gilarranz
Luis J Gilarranz 2020-1-21
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

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Map Display 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by