Setup: задайте папку данных и имя файла Служебная функция

Этот пример создает несколько временных файлов GeoTIFF и использует переменную datadir для обозначения их местоположения. Используемое здесь значение определяется выходом tempdir команда, но вы можете легко настроить это. Содержимое datadir удаляются в конце примера.

datadir = fullfile(tempdir, 'datadir');
if ~exist(datadir, 'dir')
   mkdir(datadir)
end

Задайте анонимную функцию для подготовки datadir к входу имени файла:

datafile = @(filename)fullfile(datadir, filename);

Пример 1. Запись изображения, привязанного к географическим координатам

Запись изображения, привязанного к WGS84 географическим координатам, в файл GeoTIFF. Оригинальное изображение (boston_ovr.jpg) хранится в формате JPEG со ссылкой на информацию, внешнюю по отношению к файлу изображения, в «файле привязки» (boston_ovr.jgw). Изображение обеспечивает «обзор» низкого разрешения Бостона, Массачусетса и окрестностей.

Считайте изображение из файла JPEG.

basename = 'boston_ovr';
imagefile = [basename '.jpg'];
A1 = imread(imagefile);

Получите объект привязки от файла привязки.

worldfile = getworldfilename(imagefile);
R1 = worldfileread(worldfile,'geographic',size(A1));

Запись изображения в файл GeoTIFF.

filename1 = datafile([basename '.tif']);
geotiffwrite(filename1,A1,R1)

Возврат информации о файле в виде RasterInfo объект. Обратите внимание, что значение CoordinateReferenceSystem является географической системой координат объектом.

info1 = georasterinfo(filename1);
info1.CoordinateReferenceSystem

ans = 

  geocrs with properties:

             Name: "WGS 84"
            Datum: "World Geodetic System 1984"
         Spheroid: [1×1 referenceEllipsoid]
    PrimeMeridian: 0
        AngleUnit: "degree"

Повторно импортируйте новый файл GeoTIFF и отобразите изображение обзора Бостона, правильно расположенное, в осях карты.

figure
usamap(R1.LatitudeLimits,R1.LongitudeLimits)
setm(gca,'PLabelLocation',0.05,'PlabelRound',-2,'PlineLocation',0.05)
geoshow(filename1)
title('Boston Overview')

Пример 2. Запись сетки данных, привязанной к географическим координатам

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

load topo60c
Z2 = topo60c;
R2 = topo60cR;
filename2 = datafile('topo60c.tif');
geotiffwrite(filename2,Z2,R2)

Значения в сетке данных варьируются от -7473 до 5731. Отобразите сетку как поверхность, сопоставленную с текстурой, а не как изображение интенсивности.

figure
worldmap world
gridm off
setm(gca,'MLabelParallel',-90,'MLabelLocation',90)
tmap = geoshow(filename2,'DisplayType','texturemap');
demcmap(tmap.CData)
title('Elevation Data Grid')

Пример 3: Изменение организации данных файлов GeoTIFF

Когда вы записываете данные с помощью geotiffwrite или считайте данные с помощью readgeorasterстолбцы сетки данных начинаются с севера, а строки - с запада. Например, входные данные из topo60c.mat начинается с юга, но выходные данные от topo60c.tif начинается с севера.

R2.ColumnsStartFrom
[Z3,R3] = readgeoraster(filename2);
R3.ColumnsStartFrom

ans =

    'south'


ans =

    'north'

Поэтому входные данные и данные в файле GeoTIFF перевернуты.

isequal(Z2,flipud(Z3))

ans =

  logical

   1

Если вам нужно, чтобы данные в вашей рабочей области снова совпадали, переверните значения Z и установите объект привязки так, чтобы столбцы начинались с юга:

R3.ColumnsStartFrom = 'south';
Z3 = flipud(Z3);
isequal(Z2,Z3)

ans =

  logical

   1

Данные в файле GeoTIFF закодированы с положительными значениями шкалы. Поэтому, когда вы просматриваете файл с обычным программным обеспечением просмотра TIFF, северное ребро набора данных находится в верхней части. Чтобы размещений данных в файле выхода совпадали с размещением данных входов, можно создать объект Tiff и использовать его, чтобы сбросить некоторые теги и данные изображения.

t = Tiff(filename2,'r+');

pixelScale = getTag(t,'ModelPixelScaleTag');
pixelScale(2) = -pixelScale(2);
setTag(t,'ModelPixelScaleTag',pixelScale);

tiepoint = getTag(t,'ModelTiepointTag');
tiepoint(5) = intrinsicToGeographic(R2,0.5,0.5);
setTag(t,'ModelTiepointTag',tiepoint);

setTag(t,'Compression', Tiff.Compression.None)

write(t,Z2);

rewriteDirectory(t)
close(t)

Проверьте соответствие объекта привязки и сетки данных из входных данных выходов файла данных. Для этого прочтите изображение Tiff и создайте объект ссылки. Затем сравните сетки.

t = Tiff(filename2);
Atiff = read(t);
close(t)
Rtiff = georefcells(R2.LatitudeLimits,R2.LongitudeLimits,size(Atiff));

isequal(Z2,Atiff)
isequal(R2,Rtiff)

ans =

  logical

   1


ans =

  logical

   1

Пример 4. Запись изображения, привязанного в проективную систему координат

Запишите ортофотос Concord в один файл GeoTIFF. Два смежных (с запада на восток) геореферометрических (панхроматических) ортофотоса покрывают часть Конкорда, Массачусетс, США. Файл concord_ortho.txt указывает, что данные привязаны в плановой проективной системе координат штата Массачусетс (NAD83). Модулями являются счетчики. Это соответствует коду GeoTIFF № 26986, как отмечено в спецификации GeoTIFF на http://geotiff.maptools.org/spec/geotiff6.html#6.3.3.1 под PCS_NAD83_Massachusetts.

Прочитайте два ортофота.

[X_west,R_west] = readgeoraster('concord_ortho_w.tif');
[X_east,R_east] = readgeoraster('concord_ortho_e.tif');

Объедините изображения и ссылочные объекты.

X4 = [X_west X_east];
R4 = R_west;
R4.XWorldLimits = [R_west.XWorldLimits(1) R_east.XWorldLimits(2)];
R4.RasterSize = size(X4);

Запись данных в файл GeoTIFF. Используйте номер кода 26986, указывающий на PCS_NAD83_Massachusetts Проективную систему координат.

coordRefSysCode = 26986;
filename4 = datafile('concord_ortho.tif');
geotiffwrite(filename4,X4,R4,'CoordRefSysCode',coordRefSysCode)

Возврат информации о файле в виде RasterInfo объект. Обратите внимание, что значение CoordinateReferenceSystem является проективным объектом системы координат.

info4 = georasterinfo(filename4);
info4.CoordinateReferenceSystem

ans = 

  projcrs with properties:

                    Name: "NAD83 / Massachusetts Mainland"
           GeographicCRS: [1×1 geocrs]
        ProjectionMethod: "Lambert Conic Conformal (2SP)"
              LengthUnit: "meter"
    ProjectionParameters: [1×1 map.crs.ProjectionParameters]

Отображение комбинированного ортофота Concord.

figure
mapshow(filename4)
title('Combined Orthophotos')
xlabel('MA Mainland State Plane easting, meters')
ylabel('MA Mainland State Plane northing, meters')

Пример 5: Запись обрезанного изображения из файла GeoTIFF

Запись подмножества файла GeoTIFF в новый файл GeoTIFF.

Чтение изображения RGB и ссылок на информацию из boston.tif Файл GeoTIFF.

[A5,R5] = readgeoraster('boston.tif');

Обрезать изображение.

xlimits = [ 764318  767677];
ylimits = [2951122 2954482];
[A5crop,R5crop] = mapcrop(A5,R5,xlimits,ylimits);

Запись обрезанного изображения в файл GeoTIFF. Используйте GeoKeyDirectoryTag из исходного файла GeoTIFF.

info5 = geotiffinfo('boston.tif');
filename5 = datafile('boston_subimage.tif');
geotiffwrite(filename5,A5crop,R5crop, ...
   'GeoKeyDirectoryTag',info5.GeoTIFFTags.GeoKeyDirectoryTag)

Отобразите файл GeoTIFF, содержащий обрезанное изображение.

figure
mapshow(filename5)
title('Fenway Park - Cropped Image from GeoTIFF File')
xlabel('MA Mainland State Plane easting, survey feet')
ylabel('MA Mainland State Plane northing, survey feet')

Пример 6: Запись данных о повышении в файл GeoTIFF

Запишите данные по повышению для области вокруг Южного Боулдера-Пик в Колорадо в файл GeoTIFF.

elevFilename = 'n39_w106_3arc_v2.dt1';

Считайте DEM из файла. Чтобы построить график данных с помощью geoshowданные должны иметь тип single или double. Укажите тип данных для растра, используя 'OutputType' Пара "имя-значение".

[Z6,R6] = readgeoraster(elevFilename,'OutputType','double');

Создайте структуру для хранения информации о GeoKeyDirectoryTag.

key = struct( ...
    'GTModelTypeGeoKey',[], ...
    'GTRasterTypeGeoKey',[], ...
    'GeographicTypeGeoKey',[]);

Укажите данные в системе географических координат путем определения GTModelTypeGeoKey поле как 2. Укажите, что объект ссылки использует проводки (а не камеры) путем определения GTRasterTypeGeoKey поле как 2. Укажите, что данные ссылаются на географическую систему координат путем определения GeographicTypeGeoKey поле как 4326.

key.GTModelTypeGeoKey = 2;
key.GTRasterTypeGeoKey = 2;
key.GeographicTypeGeoKey = 4326;

Запишите повышение данные в файл GeoTIFF.

filename6 = datafile('southboulder.tif');
geotiffwrite(filename6,Z6,R6,'GeoKeyDirectoryTag',key)

Проверьте, что данные записаны в файл, отобразив их. Во-первых, импортируйте векторные данные, который представляет контур состояния Колорадо с помощью shaperead. Затем отобразите контур и файл GeoTIFF.

S = shaperead('usastatelo','UseGeoCoords',true,'Selector',...
   {@(name) any(strcmp(name,{'Colorado'})),'Name'});
figure
usamap 'Colorado'
hold on
geoshow(S,'FaceColor','none')
g = geoshow(filename6,'DisplayType','mesh');
demcmap(g.ZData)
title('South Boulder Peak Elevation')

Пример 7: Запись данных без неизображений в файл TIFF

Если вы работаете с сеткой данных, которая является классом double со значениями, которые варьируются вне пределов, необходимых для изображения интенсивности с плавающей точкой (0 < = данные < = 1), и если вы храните данные в файле TIFF с помощью imwriteзатем ваши данные будут усечены до интервала [0,1], масштабированы и преобразованы в uint8. Очевидно, что некоторые или даже вся информация в исходных данных может быть потеряна. Чтобы избежать этих проблем и сохранить численный класс и область значений вашей сетки данных, используйте geotiffwrite для записи данных.

Создайте выборку данных Z.

n = 512;
Z7 = peaks(n);

Создайте объект привязки для ссылки строк и столбцов на X и Y.

R7 = maprasterref('RasterSize',[n n],'ColumnsStartFrom','north');
R7.XWorldLimits = R7.XIntrinsicLimits;
R7.YWorldLimits = R7.YIntrinsicLimits;

Создайте структуру для хранения информации о GeoKeyDirectoryTag. Установите тип модели равным 1, указывая на Проективную систему координат (PCS).

key.GTModelTypeGeoKey = 1;

Установите растровый тип равным 1, указывающий PixelIsArea (камеры).

key.GTRasterTypeGeoKey = 1;

Указать пользовательскую проективную систему координат при помощи значения 32767.

key.ProjectedCSTypeGeoKey = 32767;

Запись данных в файл GeoTIFF с geotiffwrite. Для сравнения запишите второй файл, используя imwrite.

filename_geotiff = datafile('zdata_geotiff.tif');
filename_tiff = datafile('zdata_tiff.tif');
geotiffwrite(filename_geotiff,Z7,R7,'GeoKeyDirectoryTag',key)
imwrite(Z7, filename_tiff);

Когда вы читаете файл используя imread значения данных и тип класса сохраняются. Можно увидеть, что значения данных в файле TIFF не сохранены.

geoZ  = imread(filename_geotiff);
tiffZ = imread(filename_tiff);
fprintf('Class type of Z: %s\n', class(Z7))
fprintf('Class type of data in GeoTIFF file: %s\n', class(geoZ))
fprintf('Class type of data in    TIFF file: %s\n', class(tiffZ))
fprintf('Does data in GeoTIFF file equal Z: %d\n', isequal(geoZ, Z7))
fprintf('Does data in    TIFF file equal Z: %d\n', isequal(tiffZ, Z7))

Class type of Z: double
Class type of data in GeoTIFF file: double
Class type of data in    TIFF file: uint8
Does data in GeoTIFF file equal Z: 1
Does data in    TIFF file equal Z: 0

Просмотреть сетку данных можно используя mapshow.

figure
mapshow(filename_geotiff,'DisplayType','texturemap')
title('Peaks - Stored in GeoTIFF File')

Пример 8. Изменение существующего файла при сохранении метаинформации

Можно хотеть изменить существующий файл, но сохранить большинство, если не все, метаданных в тегах TIFF. Этот пример преобразует изображение RGB из boston.tif файл в индексированное изображение и запись новых данных в индексированный файл GeoTIFF. Метаинформация TIFF, за исключением значений тегов BitDepth, BitsPerSample и PhotometricInterpretation, сохраняется.

Чтение изображения из boston.tif Файл GeoTIFF.

[A8,R8] = readgeoraster('boston.tif');

Используйте функцию MATLAB, rgb2ind, для преобразования изображения RGB в индексированное изображение X использование минимального отклонения квантования.

[X8,cmap] = rgb2ind(A8,65536);

Получите информацию о теге TIFF с помощью imfinfo.

info8 = imfinfo('boston.tif');

Создайте структуру тегов TIFF, чтобы сохранить выбранную информацию из info структура.

tags = struct( ...
    'Compression',  info8.Compression, ...
    'RowsPerStrip', info8.RowsPerStrip, ...
    'XResolution',  info8.XResolution, ...
    'YResolution',  info8.YResolution, ...
    'ImageDescription', info8.ImageDescription, ...
    'DateTime',    info8.DateTime, ...
    'Copyright',   info8.Copyright, ...
    'Orientation', info8.Orientation);

Значения для тегов PlanarConfiguration и ResolutionUnit должны быть числовыми, а не строковыми, как возвращено imfinfo. Можно задать эти теги с помощью свойств константы из класса Tiff. Например, вот возможные значения свойства константы PlanarConfiguration:

Tiff.PlanarConfiguration

ans = 

  struct with fields:

      Chunky: 1
    Separate: 2

Используйте строковое значение из info структура для получения желаемого значения.

tags.PlanarConfiguration = ...
    Tiff.PlanarConfiguration.(info8.PlanarConfiguration);

Таким же образом установите значение ResolutionUnit.

tags.ResolutionUnit = Tiff.ResolutionUnit.(info8.ResolutionUnit);

Тег Software не установлен в boston.tif файл. Однако geotiffwrite установит Software тег по умолчанию. Чтобы сохранить информацию, установите значение пустой строки, которая препятствует записи тега в файл.

tags.Software = '';

Скопируйте информацию GeoTIFF из boston.tif.

geoinfo = geotiffinfo('boston.tif');
key = geoinfo.GeoTIFFTags.GeoKeyDirectoryTag;

Запись в файл GeoTIFF.

filename8 = datafile('boston_indexed.tif');
geotiffwrite(filename8,X8,cmap,R8,'GeoKeyDirectoryTag',key,'TiffTags',tags)

Просмотр индексированного изображения.

figure
mapshow(filename8)
title('Boston - Indexed Image')
xlabel('MA Mainland State Plane easting, survey feet')
ylabel('MA Mainland State Plane northing, survey feet')

Сравните информацию в структурах, которая должна быть равной, путем печати таблицы значений.

info_rgb = imfinfo('boston.tif');
info_indexed = imfinfo(filename8);
tagNames = fieldnames(tags);
tagNames(strcmpi('Software', tagNames)) = [];
names = [{'Height' 'Width'}, tagNames'];

spacing = 2;
fieldnameLength = max(cellfun(@length, names)) + spacing;
formatSpec = ['%-' int2str(fieldnameLength) 's'];

fprintf([formatSpec formatSpec formatSpec '\n'], ...
    'Fieldname', 'RGB Information', 'Indexed Information')
fprintf([formatSpec formatSpec formatSpec '\n'], ...
    '---------', '---------------', '-------------------')

for k = 1:length(names)
    fprintf([formatSpec formatSpec formatSpec '\n'], ...
        names{k}, ...
        num2str(info_rgb.(names{k})), ...
        num2str(info_indexed.(names{k})))
end

Fieldname            RGB Information      Indexed Information  
---------            ---------------      -------------------  
Height               2881                 2881                 
Width                4481                 4481                 
Compression          Uncompressed         Uncompressed         
RowsPerStrip         256                  256                  
XResolution          300                  300                  
YResolution          300                  300                  
ImageDescription     "GeoEye"             "GeoEye"             
DateTime             2007:02:23 21:46:13  2007:02:23 21:46:13  
Copyright            "(c) GeoEye"         "(c) GeoEye"         
Orientation          1                    1                    
PlanarConfiguration  Chunky               Chunky               
ResolutionUnit       Inch                 Inch                 

Сравните информацию, которая должна быть другой, поскольку вы преобразовали изображение RGB в индексированное изображение, распечатав таблицу значений.

names = {'FileSize', 'ColorType', 'BitDepth', ...
    'BitsPerSample', 'PhotometricInterpretation'};

fieldnameLength = max(cellfun(@length, names)) + spacing;
formatSpec = ['%-' int2str(fieldnameLength) 's'];
formatSpec2 = '%-17s';

fprintf(['\n' formatSpec formatSpec2 formatSpec2 '\n'], ...
    'Fieldname', 'RGB Information', 'Indexed Information')
fprintf([formatSpec formatSpec2 formatSpec2 '\n'], ...
    '---------', '---------------', '-------------------')
for k = 1:length(names)
    fprintf([formatSpec formatSpec2 formatSpec2 '\n'], ...
        names{k}, ...
        num2str(info_rgb.(names{k})), ...
        num2str(info_indexed.(names{k})))
end

Fieldname                  RGB Information  Indexed Information
---------                  ---------------  -------------------
FileSize                   38729900         27925078         
ColorType                  truecolor        indexed          
BitDepth                   24               16               
BitsPerSample              8  8  8          16               
PhotometricInterpretation  RGB              RGB Palette      

Очистка: Удалить папку данных

Удалите временную папку и файлы данных.

rmdir(datadir, 's')

Информация о наборе данных

Файлы boston.tif и boston_ovr.jpg включает материалы, авторские права на которые принадлежат Геоглаз, все права защищены. Геоглаз была объединена в корпорацию DigitalGlobe 29 января 2013 года. Для получения дополнительной информации о наборах данных используйте команды type boston.txt и type boston_ovr.txt.

Файлы concord_orthow_w.tif и concord_ortho_e.tif получают с использованием плитки ортофото Бюро географической информации (MassGIS), Содружество Массачусетса, Административное управление по технологиям и службам безопасности. Для получения дополнительной информации о наборах данных используйте команду type concord_ortho.txt. Обновленную ссылку на данные, предоставляемые MassGIS, см. в разделе https://www.mass.gov/info-details/massgis-data-layers.

Файл DTED n39_w106_3arc_v2.dt1 предоставлено Геологической службой США.