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

Question

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

这是坐标列表：

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 功能无效。我们有足够的 RAM...

Answer 1

我建议：

% 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 可以在顶部绘制两个数据集。

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

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

Projecting lat-lon coordinates to a raster to obtain pixel values at those coordinates

matlab

raster

coordinates

map-projections