Knowledge Base

Do you have sample code to read GeoTiffs?

Sample Python Output
Sample image ouput from Python code.


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



Python

# 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

86,5 Bot