В этом примере показано, как записать данные, привязанные в стандартных географических и проективных системах координат, в файлы GeoTIFF с помощью geotiffwrite. Теговый формат графических файлов (TIFF) появился как популярный формат для хранения растровых данных. Спецификация GeoTIFF задает набор тегов TIFF, которые описывают «картографическую» информацию, связанную с растровыми данными TIFF. Используя эти теги, в файле GeoTIFF могут храниться геолокированные изображения или растровые сетки с координатами, привязанными к географической системе координат (широта и долгота) или (планарной) проективной системе координат.
Этот пример создает несколько временных файлов GeoTIFF и использует переменную datadir для обозначения их местоположения. Используемое здесь значение определяется выходом tempdir команда, но вы можете легко настроить это. Содержимое datadir удаляются в конце примера.
datadir = fullfile(tempdir, 'datadir'); if ~exist(datadir, 'dir') mkdir(datadir) end
Задайте анонимную функцию для подготовки datadir к входу имени файла:
datafile = @(filename)fullfile(datadir, filename);
Запись изображения, привязанного к 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')
Загрузка повышения растровых данных и географических камер ссылки объекта. Запись сетки данных в файл 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')
Когда вы записываете данные с помощью 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
Запишите ортофотос 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')
Запись подмножества файла 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')
Запишите данные по повышению для области вокруг Южного Боулдера-Пик в Колорадо в файл 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')
Если вы работаете с сеткой данных, которая является классом 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')
Можно хотеть изменить существующий файл, но сохранить большинство, если не все, метаданных в тегах 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 предоставлено Геологической службой США.
geotiffinfo | geotiffwrite | getworldfilename | worldfileread