In this section
  • Rendering Links ...
  • Rendering Links ...

How can I use code to read GeoTiffs?

Below the image in this article, you will find sample code in IDL, MATLAB, and Python to read in a GeoTIFF file, extract the metadata, and create an image.

The code has been tested with the following data products:

Sample Python Output

Sample image output from Python code.

IDL sample code

;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 sample code

%
% 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')

Python sample code

# 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.
#
# Required libraries:
#      numpy
#      rasterio
#      affine
#      matplotlib
#
# To use:
# Edit the file to change the GeoTIFF file name.
# Then run the code:
# python read_measures
#
import numpy as np
import rasterio
from affine import Affine
import matplotlib.pyplot as plt
import matplotlib.image as mping
import matplotlib.colors as colors

# directory path and file name
path = ''
file_in = 'greenland_vel_mosaic200_2016_2017_vel_v2.tif'

# read in file with geotiff geographic information
src = rasterio.open(path + file_in)

# print out metadata information
for k in src.meta:
  print(k,src.meta[k])

# Retrieve the affine transformation
if isinstance(src.transform, Affine):
     transform = src.transform
else:
     transform = src.affine

N = src.width
M = src.height
dx = transform.a
dy = transform.e
minx = transform.c
maxy = transform.f

# Read the image data, flip upside down if necessary
data_in = src.read(1)
if dy < 0:
  dy = -dy
  data_in = np.flip(data_in, 0)

print('Data minimum, maximum = ', np.amin(data_in), np.amax(data_in))

# Generate X and Y grid locations
xdata = minx + dx/2 + dx*np.arange(N)
ydata = maxy - dy/2 - dy*np.arange(M-1,-1,-1)

# Scale the velocities by the log of the data.
d = np.log(np.clip(data_in, 1, 3000))
data_scale = (255*(d - np.amin(d))/np.ptp(d)).astype(np.uint8)

# Construct an RGB table using a log scale between 1 and 3000 m/year.
vel = np.exp(np.linspace(np.log(1), np.log(3000), num=256))
hue = np.arange(256)/255.0
sat = np.clip(1./3 + vel/187.5, 0, 1)
value = np.zeros(256) + 0.75
hsv = np.stack((hue, sat, value), axis=1)
rgb = colors.hsv_to_rgb(hsv)
# Be sure the first color (the background) is white
rgb[0,:] = 1
cmap = colors.ListedColormap(rgb, name='velocity')

extent = [xdata[0], xdata[-1], ydata[0], ydata[-1]]
plt.figure(figsize=(8,8))
fig = plt.imshow(data_scale, extent=extent, origin='lower', cmap=cmap)
plt.title(file_in)
# Hide the axes and remove the space around them
plt.axis('off')
fig.axes.get_xaxis().set_visible(False)
fig.axes.get_yaxis().set_visible(False)
tickval = [1, 10, 100, 1000, 3000]
t = np.log(tickval)
cb = plt.colorbar(fig, ticks=255*(t - t[0])/(t[-1] - t[0]), shrink=0.5)
cb.set_label('Velocity Magnitude (m/year)')
cb.ax.set_yticklabels(tickval)
plt.savefig("python_sample.png", dpi=300, bbox_inches='tight', pad_inches=0.5)
plt.show()
 

Last Updated June 06, 2018

 

 

Related Data Collections

MEaSUREs