-
Notifications
You must be signed in to change notification settings - Fork 0
/
SumArea.m
65 lines (47 loc) · 3.05 KB
/
SumArea.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
function [TotalArea,FreeArea,Lats,Lons] = SumArea(grid, gridR, combW)
%SumArea calculates the total area for two cases
% Calculates the total area for the two following cases:
% 1. The total area bounded within the borders of the contiguous US
% 2. The land area remaining after the following map layers have been
% subtracted from the full US map:
% - agricultural development
% - standing water
% - urban development
% - forest coverage
% - elevation gradient
fprintf('\n------------------------') % Print current function to screen
fprintf('\n---Area Calculations----\n')
LL = fprintf('Complete: XXX.XXX');
latvec = flip(linspace(gridR.LatitudeLimits(1), ... % Generate vector of reference latitudes
gridR.LatitudeLimits(2),size(grid,1)));
lonvec = linspace(gridR.LongitudeLimits(1),gridR.LongitudeLimits(2), ...% Generate vector of reference longitudes
size(grid,2));
earthellipsoid = referenceSphere('earth','mile'); % Set the reference datum for the area calculations
TotalArea = 0; % Initialize total US area variable
FreeArea = 0; % Initialize total wind energy development area variable
Total = length(latvec) * length(lonvec); % Total number of quadrants in bounded area
Count = 0; % Progress counter
Tally = 1; % Free quadrant counter
for j = 2:size(grid,2)
for i = 2:size(grid,1)
if mod(Count,10000) == 0
fprintf(repmat('\b',1,LL))
fprintf('Complete: %6.2f%%',(Count/Total)*100) % Print progress to screen
end
if ~isnan(grid(i,j))
TotalArea = TotalArea + areaquad(latvec(i-1), ... % Sum total area in the continental US
lonvec(j-1),latvec(i),lonvec(j),earthellipsoid);
end
if ~isnan(combW(i,j))
FreeArea = FreeArea + areaquad(latvec(i-1),lonvec(j-1), ... % Sum area available for wind energy expansion
latvec(i),lonvec(j),earthellipsoid);
Lats(Tally) = i; % Record latitudes of available wind farm site
Lons(Tally) = j; % Record longitudes of available wind farm site
Tally = Tally + 1; % Tally number of sites
end
Count = Count + 1; % Update progress counter
end
end
fprintf(repmat('\b',1,LL))
fprintf('------------------------\n') % Print progress to screen
end