MATLAB Answers

USGS global landcover in lat/long format?

4 views (last 30 days)
Bryan
Bryan on 2 May 2014
Edited: Bryan on 6 May 2014
Hi,
I would like to load the lat/long Global Landcover characterisation file from USGS.
MatLab provides tools for loading the USGS landcover data in Goode projection (function avhrrgoode) and Lambert projection (function avhrrlambert) but not lat/long projection.
The file simply is binary file representing a large raster of the Earth using a 30 arcsecond grid. The file is described under section 9.2 in the following page: edc2.usgs.gov/glcc/globdoc2_0.php
I have tried opening using fread (runs out of memory) and multibandread (freezes the computer) but so far I have had no joy.
Anybody have experience with this?
Thanks Bryan

  2 Comments

dpb
dpb on 2 May 2014
No, no experience but look at
doc memmapfile

Sign in to comment.

Accepted Answer

Bryan
Bryan on 5 May 2014
Edited: Bryan on 6 May 2014
Thanks very much to dpb for making me aware of the MatLab function memmapfile, it set me on the right path.
I wrote a script for opening the file in case anybody ends up running into the same problem in future and finds this discussion by searching the MatLab questions. I cleaned it up a bit, see below:
% Crop USGS Global Landcover lat/lon file in half-minute resolution
% From Global Land Cover Characteristics Data Base Version 2.0
% MatLab script by Bryan Lougheed
% Cropped window returned as matrix "landcover"
% Specify location of file
file_loc = ['C:\GISdata\gusgs2_0ll.img'];
% Specify window to be cropped below
% In arc seconds, must coincide exactly with USGS arc-second coordinates
% values are the coordinates the col/row is centred on
% window in arc seconds
west_lon_as = -10.5 *60*60+15; % west column lon (minus for W, plus for E)
east_lon_as = -5.4 *60*60+15; % east column lon (minus for W, plus for E)
north_lat_as = +56 *60*60+15; % north row lat (minus for S, plus for N)
south_lat_as = +51 *60*60+15; % south row lat (minus forS, plus for N)
%---now get data-----------------------------------------------------------
%load file & display file info in command window
usgs = memmapfile(file_loc)
% Arc second coordinates in USGS file
% See Section 9.2 here: http://edc2.usgs.gov/glcc/globdoc2_0.php
all_lons_as = -647985:30:647985;
all_lats_as = 323985:-30:-323985;
% col indexes
start_col = find(all_lons_as == west_lon_as);
end_col = find(all_lons_as == east_lon_as);
% row indexes
start_row = find(all_lats_as == north_lat_as);
end_row = find(all_lats_as == south_lat_as);
[X,Y] = meshgrid( west_lon_as:30:east_lon_as , north_lat_as:-30:south_lat_as );
%load from file line by line
newrow = 0;
landcover = NaN(size(X));
for i = start_row:1:end_row;
start_index = ((i-1)*360*60*2)+start_col;
end_index = ((i-1)*360*60*2)+end_col;
newrow = newrow+1;
landcover(newrow,:) = usgs.Data(start_index:end_index);
end
% preview the window
image(landcover)
% convert meshgrid to decimal degrees
X = X ./ 3600;
Y = Y ./ 3600;

  0 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by