-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathlandfind_indices.m
148 lines (98 loc) · 5.63 KB
/
landfind_indices.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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
function [landmask_ind] = landfind_indices(longwest,longeast,latsouth,latnorth,pixel_dim)
% function that outputs land mask indices, given four coordinate bounds (MUST be the same as the min and max of the image being overlaid) and the pixel dimensions
%begin bathymetry calculations
ncid_topo = netcdf.open('topo2.grd','NC_NOWRITE');
% x_range = netcdf.getVar(ncid_topo,0);
% y_range = netcdf.getVar(ncid_topo,1);
% z_range = netcdf.getVar(ncid_topo,2);
% spacing = netcdf.getVar(ncid_topo,3);
% dimension = netcdf.getVar(ncid_topo,4);
z = netcdf.getVar(ncid_topo,5);
spacing = 1/30; %spacing of bathymetry data, in degrees
%masking bathymetry threshold (in meters elevation)
mask_threshold = 0;
% %isolate bathymetry data
%
% bathy_ns_indices = (floor((90-latnorth)/spacing)+1:ceil((90-latsouth)/spacing))';
%
% lat_bathy = 90 - (double(bathy_ns_indices-1)*spacing + spacing/2);
%
% bathy_ew_indices = (floor(longwest/spacing)+1:ceil(longeast/spacing))';
%
% long_bathy = double(bathy_ew_indices-1)*spacing + spacing/2;
lat_bathy = 90 - ((0.5*spacing):spacing:180)';
long_bathy = ((0.5*spacing):spacing:360)';
z_reshaped = reshape(z,360/spacing,180/spacing);
clear z
% z_range = z_reshaped(bathy_ew_indices,bathy_ns_indices);
%
% clear z_reshaped
%
% z_range = flipdim(z_range,2);
z_reshaped = flipdim(z_reshaped,2);
lat_bathy = flipdim(lat_bathy,1);
mask_ind = find(z_reshaped >= mask_threshold);
zero_padding_dim = [5 5];
z_land_mask = zeros(size(z_reshaped));
% mask_ind_padded = mask_ind + (zero_padding_dim(2)*(size(z_reshaped,1) + (2*zero_padding_dim(1)))) + zero_padding(1);
z_land_mask(mask_ind) = 1;
% resolutions of imported image
% x_res_image = (longeast - longwest)/(pixel_dim(1) - 1);
% y_res_image = (latnorth - latsouth)/(pixel_dim(2) - 1);
x_res_image = (mod(longeast - longwest - 1e-10,360) + 1e-10)/(pixel_dim(1));
y_res_image = (latnorth - latsouth)/(pixel_dim(2));
x_res_ratio = x_res_image/spacing;
y_res_ratio = y_res_image/spacing;
% apply smoothing filter
x_smooth_filter_vec = zeros(size(z_land_mask,1),1);
y_smooth_filter_vec = zeros(size(z_land_mask,2),1);
edge_full_ind = floor((x_res_ratio + 1)/2);
x_smooth_filter_vec(1:edge_full_ind) = 1;
x_smooth_filter_vec(edge_full_ind + 1) = ((x_res_ratio + 1)/2) - edge_full_ind;
x_smooth_filter_vec((length(x_smooth_filter_vec) - edge_full_ind + 1):length(x_smooth_filter_vec)) = flipdim(x_smooth_filter_vec(2:(edge_full_ind + 1)),1);
y_smooth_filter_vec(1:edge_full_ind) = 1;
y_smooth_filter_vec(edge_full_ind + 1) = ((y_res_ratio + 1)/2) - edge_full_ind;
y_smooth_filter_vec((length(y_smooth_filter_vec) - edge_full_ind + 1):length(y_smooth_filter_vec)) = flipdim(y_smooth_filter_vec(2:(edge_full_ind + 1)),1);
x_smooth_filter_array = repmat(x_smooth_filter_vec,[1 size(z_land_mask,2)]);
y_smooth_filter_array = repmat(y_smooth_filter_vec',[size(z_land_mask,1) 1]);
smooth_filter_array = x_smooth_filter_array.*y_smooth_filter_array;
smooth_filter_array = smooth_filter_array/(sum(sum(smooth_filter_array)));
fft_smooth_filter_array = fft(fft(smooth_filter_array,[],1),[],2);
z_land_mask_smoothed = ifft(ifft(fft_smooth_filter_array.*fft(fft(z_land_mask,[],1),[],2),[],2),[],1);
% "interpolate" to the same grid as the image
if mod(longwest,360) >= mod(longeast,360)
in_bathy_lon_range_ind_1 = find(long_bathy >= mod(longwest,360));
in_bathy_lon_range_ind_2 = find(long_bathy <= mod(longeast,360));
in_bathy_lon_range_ind = [in_bathy_lon_range_ind_1; in_bathy_lon_range_ind_2];
else
in_bathy_lon_range_ind = find((long_bathy >= mod(longwest,360)) & (long_bathy <= mod(longeast,360)));
end
in_bathy_lat_range_ind = find((lat_bathy >= latsouth) & (lat_bathy <= latnorth));
% SW_corner_x_weight_inner = (longwest - long_bathy(in_bathy_lon_range_ind(1) - 1))/spacing;
% SW_corner_y_weight_inner = (latsouth - lat_bathy(min(in_bathy_lat_range_ind) - 1))/spacing;
SW_corner_x_weight_inner = ((mod(longwest - long_bathy(in_bathy_lon_range_ind(1)) + 180,360) - 180)/spacing) + 1;
SW_corner_y_weight_inner = ((latsouth - lat_bathy(min(in_bathy_lat_range_ind)))/spacing) + 1;
x_rel_placing_vec = SW_corner_x_weight_inner + (x_res_ratio*((0:1:(pixel_dim(1) - 1))'));
land_mask_x_interp = NaN(pixel_dim(1),size(z_land_mask_smoothed,2));
for x_ind = 1:length(x_rel_placing_vec)
curr_x_rel_placing = x_rel_placing_vec(x_ind);
x_bathy_lower_ind = mod(in_bathy_lon_range_ind(1) - 1 + floor(curr_x_rel_placing) - 1,length(long_bathy)) + 1;
x_bathy_upper_ind = mod(x_bathy_lower_ind + 1 - 1,length(long_bathy)) + 1;
x_upper_weight = curr_x_rel_placing - floor(curr_x_rel_placing);
x_lower_weight = 1 - x_upper_weight;
land_mask_x_interp(x_ind,:) = (x_lower_weight*z_land_mask_smoothed(x_bathy_lower_ind,:)) + (x_upper_weight*z_land_mask_smoothed(x_bathy_upper_ind,:));
end
y_rel_placing_vec = SW_corner_y_weight_inner + (y_res_ratio*((0:1:(pixel_dim(2) - 1))'));
land_mask = NaN(pixel_dim);
for y_ind = 1:length(y_rel_placing_vec)
curr_y_rel_placing = y_rel_placing_vec(y_ind);
y_bathy_lower_ind = max([1 (min(in_bathy_lat_range_ind) - 1 + floor(curr_y_rel_placing))]);
y_bathy_upper_ind = min([length(lat_bathy) (y_bathy_lower_ind + 1)]);
y_upper_weight = curr_y_rel_placing - floor(curr_y_rel_placing);
y_lower_weight = 1 - y_upper_weight;
land_mask(:,y_ind) = (y_lower_weight*land_mask_x_interp(:,y_bathy_lower_ind)) + (y_upper_weight*land_mask_x_interp(:,y_bathy_upper_ind));
end
% find indices of land points
landmask_ind = find(land_mask > 0.5);
% close netCDF file
netcdf.close(ncid_topo);