Проецирование координат широты на растр для получения значений пикселей в этих координатах - PullRequest
Новая
       74

Проецирование координат широты на растр для получения значений пикселей в этих координатах

0 голосов
/ 21 ноября 2018

У нас есть растровый файл с географической привязкой и список координат широты и долготы, которые мы хотим построить поверх этого растрового файла, чтобы получить растровые значения по этим координатам.

Вот список координат: 

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

enter image description here

Конечно, долгота находится между -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)

enter image description here

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

В результате на этом рисунке enter image description here

Проблема в том, что проекционные координаты идут от -3 до +3 и от -1 до +1и ось в растровом файле изменяется от -2 до +2, но до степени 7, поэтому, если мы нанесем эти точки поверх этого растра, мы увидим все как одну единственную точку где-нибудь в Атлантике.

В конце концов, мы хотим использовать функцию latlon2pix, чтобы получить значения пикселей в каждом из tон координирует точки, но, во-первых, мы должны быть способны соединить две вещи.Есть идеи?

Функция геошоу не работает.У нас достаточно оперативной памяти ...

1 Ответ

0 голосов
/ 05 декабря 2018

Я предлагаю следующее: 

% 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');
...