时间戳

首页 > 解决方案 > 将经纬度坐标投影到栅格以获得这些坐标处的像素值

matlab - 将经纬度坐标投影到栅格以获得这些坐标处的像素值

问题描述

我们有一个地理参考的栅格文件和一个纬度-经度坐标列表，我们希望在该栅格文件的顶部绘制这些坐标以获得这些坐标处的栅格值。

这是坐标列表：

mapshow(lon,lat,'DisplayType','point')  %x,y

在此处输入图像描述

当然，经度介于 -180 和 +180 之间，纬度介于 -90 和 +90 之间。

这是光栅文件：

proj = geotiffinfo('global2001.tif');
[raster_tiff,cmap_tiff,reference_tiff] = geotiffread('global2001.tif');
figure
mapshow(raster_tiff,cmap_tiff,reference_tiff)

在此处输入图像描述

proj =           FileModDate: '20-nov-2018 14:15:52'
             FileSize: 121625752
               Format: 'tif'
        FormatVersion: []
               Height: 40032
                Width: 80062
             BitDepth: 8
            ColorType: 'indexed'
            ModelType: 'ModelTypeProjected'
                  PCS: ''
           Projection: ''
               MapSys: ''
                 Zone: []
         CTProjection: 'CT_Sinusoidal'
             ProjParm: [7×1 double]
           ProjParmId: {7×1 cell}
                  GCS: 'WGS 84'
                Datum: 'World Geodetic System 1984'
            Ellipsoid: 'WGS 84'
            SemiMajor: 6378137
            SemiMinor: 6.3568e+06
                   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]

现在我们使用与栅格文件相同的投影来投影坐标：

mstruct = geotiff2mstruct(proj);
% get current axis 
axesm('MapProjection',mstruct.mapprojection,'Frame','on')
h=get(gcf,'CurrentAxes');
assert(ismap(h)==1,'Current axes must be map axes.')
mstruct=gcm;
[x,y] = mfwdtran(mstruct,lat,lon,h,'none');

figure
mapshow(x,y,'DisplayType','point')  %x,y

导致这个图 在此处输入图像描述

问题是投影坐标从-3到+3，从-1到+1，光栅文件中的轴从-2到+2但是是7的幂，所以如果我们把这些点画在上面在该栅格中，我们将所有内容视为大西洋某处的一个点。

最终，我们想使用函数 latlon2pix 来获取每个坐标点的像素值，但首先，我们需要能够将这两个东西放在一起。有任何想法吗？

geoshow 功能不起作用。我们有足够的内存...

标签： matlabcoordinatesrastermap-projections

解决方案

我建议这样做：

% assigning reference for tif file    
proj = geotiffinfo('global2001.tif');
% reading the tif into raster_tiff, storing the colormap, save the meta data 
[raster_tiff,cmap_tiff,reference_tiff] = geotiffread('global2001.tif');

投影可以使用projfwd完成，投影矩阵存储在变量 proj 中。然后，数据点将以栅格文件的（此处为正弦）地图坐标表示。

% using the projection matrix (proj.RefMatrix) and computing the raster coordinates in meters
[x,y] = projfwd(proj,lat,lon);

使用mapshow one 可以在顶部绘制两个数据集。

mapshow(raster_tiff,cmap_tiff,reference_tiff);
mapshow(x,y,'DisplayType', 'point', 'Color', 'm',...
            'MarkerEdgeColor', 'auto');

推荐阅读