USGS global landcover in lat/long format?

8 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

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;

More Answers (0)

Categories

Find more on Weather and Atmospheric Science in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!