Настройка: задайте служебную функцию папки и имени файла данных

Этот пример создает несколько временных файлов 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)

Заметьте, что информация о GeoTIFF из файла указывает, что Географическая система координат (GCS) является "WGS 84" и что все поля (PCS, Проекция, MapSys, Зона, CTProjection, ProjParm, ProjParmId и UOMLength) сопоставленный со спроектированной системой координат пусты.

geotiffinfo(filename1)

ans = 

  struct with fields:

             Filename: 'C:\Users\jbenham\AppData\Local\Temp\datadir\boston_ovr.tif'
          FileModDate: '06-Jan-2020 09:40:21'
             FileSize: 1674212
               Format: 'tif'
        FormatVersion: []
               Height: 769
                Width: 722
             BitDepth: 8
            ColorType: 'truecolor'
            ModelType: 'ModelTypeGeographic'
                  PCS: ''
           Projection: ''
               MapSys: ''
                 Zone: []
         CTProjection: ''
             ProjParm: []
           ProjParmId: ''
                  GCS: 'WGS 84'
                Datum: 'World Geodetic System 1984'
            Ellipsoid: 'WGS 84'
            SemiMajor: 6378137
            SemiMinor: 6.3568e+06
                   PM: 'Greenwich'
    PMLongToGreenwich: 0
            UOMLength: ''
    UOMLengthInMeters: []
             UOMAngle: 'degree'
    UOMAngleInDegrees: 1
            TiePoints: [1×1 struct]
           PixelScale: [3×1 double]
           SpatialRef: [1×1 map.rasterref.GeographicCellsReference]
            RefMatrix: [3×2 double]
          BoundingBox: [2×2 double]
         CornerCoords: [1×1 struct]
         GeoTIFFCodes: [1×1 struct]
          GeoTIFFTags: [1×1 struct]

Повторно импортируйте новый файл 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 topo
Z2 = topo;
R2 = georefcells(topolatlim,topolonlim,size(Z2));
filename2 = datafile('topo.tif');
geotiffwrite(filename2,Z2,R2)

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

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

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

Когда вы записываете данные с помощью geotiffwrite или считайте данные с помощью readgeoraster, столбцы сетки данных запускаются с севера, и строки запускаются с запада. Например, входные данные от topo.mat запускается с юга, но выходных данных от topo.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: запишите изображение, ссылаемое в спроектированную систему координат

Запишите ортофотографии Согласия в один файл 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)

Заметьте, что информация о GeoTIFF из файла указывает, что Географическая система координат (GCS) является "NAD83" и что все поля, сопоставленные со Спроектированной Системой координат, заполнены в соответствующими значениями.

geotiffinfo(filename4)

ans = 

  struct with fields:

             Filename: 'C:\Users\jbenham\AppData\Local\Temp\datadir\concord_ortho.tif'
          FileModDate: '06-Jan-2020 09:40:26'
             FileSize: 8068660
               Format: 'tif'
        FormatVersion: []
               Height: 2000
                Width: 4000
             BitDepth: 8
            ColorType: 'grayscale'
            ModelType: 'ModelTypeProjected'
                  PCS: 'NAD83 / Massachusetts Mainland'
           Projection: 'SPCS83 Massachusetts Mainland zone (meters)'
               MapSys: 'STATE_PLANE_83'
                 Zone: 2001
         CTProjection: 'CT_LambertConfConic_2SP'
             ProjParm: [7×1 double]
           ProjParmId: {7×1 cell}
                  GCS: 'NAD83'
                Datum: 'North American Datum 1983'
            Ellipsoid: 'GRS 1980'
            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]

Отобразите объединенные ортофотографии Согласия.

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 поле как 4 326.

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

Если вы работаете с сеткой данных, которая является классом дважды со значениями, которые располагаются вне пределов, требуемых изображения интенсивности с плавающей точкой (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;

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

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);

Тег программного обеспечения не установлен в 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 включайте материалы, защищенные авторским правом GeoEye, все права защищены. GeoEye был объединен в корпорацию DigitalGlobe 29-го января 2013. Для получения дополнительной информации о наборах данных, используйте команды type boston.txt и type boston_ovr.txt.

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

Файл DTED n39_w106_3arc_v2.dt1 любезность Геологической службы США.