USGS global landcover in lat/long format?

3 次查看（过去 30 天）
Bryan
Bryan 2014-5-2
编辑： Bryan 2014-5-6
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

请先登录，再回答此问题。

采纳的回答

Bryan
Bryan 2014-5-5
编辑：Bryan 2014-5-6
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 个）

类别

在 Help Center 和 File Exchange 中查找有关 Data Import and Analysis 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by