В этом примере показано, как записывать данные, связанные со стандартными географическими и проекционными системами координат, в файлы GeoTIFF с помощью geotiffwrite. Формат Tagged-Image File Format (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 Спроецированную систему координат (Projected Coordinate System).
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. Используйте тег GeoKeyStartTag из исходного файла 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');
Создайте структуру, в которой будет храниться информация о GeoKeyStartTag.
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 < = data < = 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;
Создайте структуру, в которой будет храниться информация о GeoKeyStartTag. Задайте тип модели 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 и ResolityUnit должны быть числовыми, а не строковыми, как возвращено imfinfo. Эти теги можно задать с помощью свойств константы класса Tiff. Например, ниже приведены возможные значения для свойства константы PlanarConfiguration:
Tiff.PlanarConfiguration
ans =
struct with fields:
Chunky: 1
Separate: 2
Используйте строковое значение из info структура для получения требуемого значения.
tags.PlanarConfiguration = ...
Tiff.PlanarConfiguration.(info8.PlanarConfiguration);
Установите значение ResolityUnit таким же образом.
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 включают материалы, авторские права на которые принадлежат GeoEye, все права защищены. GeoEye была объединена в корпорацию DigityGlobe 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