Knowledge Base

Do you have sample code to read GeoTiffs?

Below you will find  sample code in IDL and MatLab to read in a GeoTIFF file, extract the metadata, and create an image.
The code has been tested with data from:
MEaSUREs Greenland Ice Sheet Velocity Mosaics SAR and Landsat, Version 1 (NSIDC-0725) https://doi.org/10.5067/OBXCG75U7540
and
MEaSUREs Greenland Ice Sheet Velocity Map InSAR Data, Version 2 (NSIDC-0478) https://doi.org/10.5067/OC7B04ZM9G6Q

IDL

;pro read_geotiff
;
; NSIDC, April 2018.
;
; Sample program to read in a GeoTIFF file, extract the metadata,
; and create an image.
;
; The code has been tested with data from:
;   NSIDC-0725, MEaSUREs Greenland Ice Sheet Velocity Mosaics SAR and Landsat, Version 1
;     https://doi.org/10.5067/OBXCG75U7540
;   NSIDC-0478, MEaSUREs Greenland Ice Sheet Velocity Map InSAR Data, Version 2
;     https://doi.org/10.5067/OC7B04ZM9G6Q
;
; Note: Custom color map designed by Terry Haran, NSIDC.
;
; To use:
; Edit the file to change the GeoTIFF file name.
; Then run the code:
; IDL> .run read_measures
;
; directory path and file name
path = FILE_DIRNAME(ROUTINE_FILEPATH(), /MARK)
file_in = 'greenland_vel_mosaic500_2016_2017_vel_v2.tif'

; read in file with geotiff geographic information
data_in = READ_TIFF(path + file_in, GEOTIFF=GeoKeys, ORIENTATION=orient)
if (orient eq 1) then begin
  data_in = REVERSE(data_in, 2, /OVERWRITE)
endif

; print out metadata information
array_info = SIZE(data_in, /STRUCTURE)
print, 'IDL datatype = ', array_info.TYPE_NAME
print, 'Number of elements = ', array_info.N_ELEMENTS
print, 'Number of dimensions = ', array_info.N_DIMENSIONS
print, 'Dimensions = ', array_info.DIMENSIONS[0:array_info.N_DIMENSIONS-1]
minn = MIN(data_in, MAX=maxx)
print, 'Data minimum, maximum = ', minn, maxx
foreach value, Hash(GeoKeys), key do begin
  print, key, '...', value
endforeach

N     = array_info.DIMENSIONS[0]
M     = array_info.DIMENSIONS[1]
dx    = GeoKeys.ModelPixelScaleTag[0]
dy    = GeoKeys.ModelPixelScaleTag[1]
minx  = GeoKeys.ModelTiepointTag[3]
maxy  = GeoKeys.ModelTiepointTag[4]

; Generate X and Y grid locations
xdata = minx + dx/2 + ([0:N-1]*dx)
ydata = maxy - dy/2 - ([M-1:0:-1]*dy)

; Scale the velocities by the log of the data.
data_scale = BYTSCL(alog(data_in > 1 < 3000))

; Construct an RGB table using a log scale between 1 and 3000 m/year.
vel = exp(INTERPOL([alog(1), alog(3000)], 256))
hue = [0:255]/255.0*360
sat = (1./3 + vel/187.5) > 0 < 1
value = REPLICATE(0.75, 256)
COLOR_CONVERT, hue, sat, value, r, g, b, /HSV_RGB
ctable = [[r], [g], [b]]
; Be sure the first color (the background) is white
ctable[0,*] = 255

i = IMAGE(data_scale, xdata, ydata, DIM=[800,800], $
  RGB_TABLE=ctable, TITLE=file_in, $
  GEOTIFF=GeoKeys, LIMIT=[60, -80,85,-10])
grid = MAPGRID(GRID_LONGITUDE=10, GRID_LATITUDE=5)
; MAPGRID will extend the map out to the grid boundary,
; so reset the x/yrange to clip the grid.
i.xrange = [min(xdata), max(xdata)]
i.yrange = [min(ydata), max(ydata)]
tickval = [1, 10, 100, 1000, 3000]
c = COLORBAR(TARGET=i, ORIENTATION=1, TEXTPOS=1, $
  TITLE='Velocity Magnitude (m/year)', SUBTICKLEN=0, $
  TICKVALUES=BYTSCL(ALOG(tickval)), TICKNAME=STRTRIM(tickval,2))

i.Save, 'idl_sample.png', RESOLUTION=300, BORDER=20

end

Matlab

%
% NSIDC, April 2018.
%
% Sample program to read in a GeoTIFF file, extract the metadata,
% and create an image.
%
% The code has been tested with data from:
%   NSIDC-0725, MEaSUREs Greenland Ice Sheet Velocity Mosaics SAR and Landsat, Version 1
%     https://doi.org/10.5067/OBXCG75U7540
%   NSIDC-0478, MEaSUREs Greenland Ice Sheet Velocity Map InSAR Data, Version 2
%     https://doi.org/10.5067/OC7B04ZM9G6Q
%
% Note: Custom color map designed by Terry Haran, NSIDC.
%
% Edit the file to change the GeoTIFF file name.
% Then run the code within the Matlab application.
%
% directory path and file name
path = '';
file_in = 'greenland_vel_mosaic500_2016_2017_vel_v2.tif';

% Retrieve the image info
GeoKeys = imfinfo(strcat(path,file_in))
N     = GeoKeys.Width;
M     = GeoKeys.Height;
dx    = GeoKeys.ModelPixelScaleTag(1);
dy    = GeoKeys.ModelPixelScaleTag(2);
minx  = GeoKeys.ModelTiepointTag(4);
maxy  = GeoKeys.ModelTiepointTag(5);

% Generate x/y pixel location vectors
xdata = minx + dx/2 + ((0:N-1).*dx);
ydata = maxy - dy/2 - ((M  -1:-1:0).*dy);

% Read the image data, flip upside down if necessary
data_in = imread(strcat(path,file_in));
if GeoKeys.Orientation
    data_in = flipud(data_in);
end

% Clip the data and scale the log of the data to a byte range.
data_in(data_in > 3000) = 3000;
data_in(data_in < 1) = 1;
d = log(data_in);
dmin = min(d(:));
dmax = max(d(:));
data_scale = (d - dmin)./(dmax - dmin).*255;
data_scale = uint8(data_scale);

% Construct an RGB table using a log scale between 1 and 3000 m/year.
vel = exp(linspace(log(1), log(3000), 256));
hue = (0:255)/255.0;
sat = 1./3 + vel/187.5;
sat(sat < 0) = 0;
sat(sat > 1) = 1;
value = 0.75*ones(1,256);
ctable = hsv2rgb(cat(2, hue', sat', value'));
% Be sure the first color (the background) is white
ctable(1,:) = 1;

figure('pos', [10,10,800,800])
i = imagesc(xdata,ydata,data_scale);
title(file_in, 'Interpreter', 'none', 'FontSize', 14)
axis equal xy; axis image off
tickval = [1, 10, 100, 1000, 3000];
t = log(tickval);
tmin = min(t(:));
tmax = max(t(:));
t = (t - tmin)./(tmax - tmin).*255;
h = colorbar('Ticks', t,...
         'TickLabels', tickval, 'FontSize', 12);
set(get(h,'label'), 'string', 'Velocity Magnitude (m/year)');
set(h,'position',[.8 .35 .05 .3]);
colormap(ctable)

saveas(i, 'matlab_sample.png')
 

Last Updated May 07, 2018

86,5 Bot