USGS global landcover in lat/long format?

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

[dpb]

No, no experience but look at

doc memmapfile

[Bryan]

Thanks!

Answer 1

[Bryan]

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;