-
Notifications
You must be signed in to change notification settings - Fork 97
/
EarthView.m
314 lines (270 loc) · 11.4 KB
/
EarthView.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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
%EarthView Image from Google maps
%
% A concrete subclass of ImageSource that acquires images from Google
% maps.
%
% Methods::
% grab Grab a frame from Google maps
% size Size of image
% close Close the image source
% char Convert the object parameters to human readable string
%
% Examples::
% Create an EarthView camera
% ev = EarthView();
% Zoom into QUT campus in Brisbane
% ev.grab(-27.475722,153.0285, 17);
% Show aerial view of Brisbane in satellite and map view
% ev.grab('brisbane', 14)
% ev.grab('brisbane', 14, 'map')
%
% Notes::
% - Google limit the number of map queries limit to 1000 unique (different)
% image requests per viewer per day. A 403 error is returned if the daily
% quota is exceeded.
% - Maximum size is 640x640 for free access, business users can get more.
% - There are lots of conditions on what you can do with the images,
% particularly with respect to publication. See the Google web site for
% details.
%
% Author::
% Peter Corke, with some lines of code from from get_google_map by Val
% Schmidt.
%
% See also ImageSource.
%TODO:
% map features/element control
% street view
% roadmap, satellite, terrain, hybrid
classdef EarthView < ImageSource
% e = EarthView() is an object that returns images of the Earth's surface
% obtained from Google Maps.
%
% e = EarthView(googlekey) as above but the google key string is passed in.
%
% e.grab(lat, lon, zoom, options)
properties
key
type
scale
end
methods
function ev = EarthView(varargin)
%EarthView.EarthView Create EarthView object
%
% EV = EarthView(OPTIONS)
%
% Options::
% 'satellite' Retrieve satellite image
% 'map' Retrieve map image
% 'hybrid' Retrieve satellite image with map overlay
% 'scale' Google map scale (default 18)
% 'width',W Set image width to W (default 640)
% 'height',H Set image height to H (default 640)
% 'key',S The Google maps key string
%
% see also options for ImageSource.
%
% Notes::
% - A key is required before you can use the Google Static Maps API.
% The key is a long string that can be passed to the constructor or
% saved as an environment variable GOOGLE_KEY. You need a Google
% account before you can register for a key.
%
% Notes::
% - Scale is 1 for the whole world, ~20 is about as high a resolution
% as you can get.
%
% See also ImageSource, EarthView.grab.
%TODO
% clone this for new StreetView API
% method to return lat/long or NE matrices corresp to pixels
ev = ev@ImageSource(varargin{:});
opt.type = {'satellite', 'map', 'hybrid'};
opt.scale = 18;
opt.key = [];
[opt,args] = tb_optparse(opt, varargin);
ev.type = opt.type;
ev.scale = 1;
% set default size params if not set
if isempty(ev.width)
ev.width = 640;
ev.height = 640;
end
if isempty(opt.key)
opt.key = getenv('GOOGLE_KEY');
end
ev.key = opt.key;
end
function [im,E,N] = grab(ev, varargin)
% EarthView.grab Grab an aerial image
%
% IM = EV.grab(LAT, LONG, OPTIONS) is an image of the Earth
% centred at the geographic coordinate (lat, long).
%
% IM = EarthView.grab(LAT, LONG, ZOOM, OPTIONS) as above with the
% specified zoom. ZOOM is an integer between 1 (zoom right out)
% to a maximum of 18-20 depending on where in the world you
% are looking.
%
% [IM,E,N] = EarthView.grab(LAT, LONG, OPTIONS) as above but also
% returns the estimated easting E and northing N. E and N are both
% matrices, the same size as IM, whose corresponding elements are
% the easting and northing are the coordinates of the pixel.
%
% [IM,E,N] = EarthView.grab(NAME, OPTIONS) as above but uses a
% geocoding web site to resolve the name to a location.
%
% Options::
% 'satellite' Retrieve satellite image
% 'map' Retrieve map image
% 'roadmap' Retrieve map image (synonym for 'map')
% 'terrain' Retrieve terrain map
% 'hybrid' Retrieve satellite image with map overlay
% 'roads' Retrieve a binary image that shows only roads, no labels.
% Roads are white, everything else is black.
% 'noplacenames' Don't show placenames.
% 'noroadnames' Don't show roadnames.
%
%
% Examples::
% Zoom into QUT campus in Brisbane
% ev.grab(-27.475722,153.0285, 17);
% Show aerial view of Brisbane in satellite and map view
% ev.grab('brisbane', 14)
% ev.grab('brisbane', 14, 'map')
%
% Notes::
% - If northing/easting outputs are requested the function
% deg2utm is required (from MATLAB Central)
% - The easting/northing is somewhat approximate, see
% get_google_map on MATLAB Central.
% - If no output argument is given the image is displayed using idisp.
opt.type = {'satellite', 'map', 'hybrid', 'terrain', 'roadmap', 'roads'};
opt.scale = ev.scale;
opt.roadnames = true;
opt.placenames = true;
opt.landscape = false;
opt.onlyroads = false;
[opt,args] = tb_optparse(opt, varargin);
switch opt.type
case 'map'
opt.type = 'roadmap'
case 'roads'
opt.type = 'roadmap';
opt.onlyroads = true;
end
% build the URL
if ischar(args{1})
% given a string name, do a geocode lookup
place = args{1};
zoom = args{2};
% build the URL, and load the XML document
url = sprintf('http://maps.googleapis.com/maps/api/geocode/xml?address=%s&sensor=false', place);
doc = xmlread(url);
% walk the XML document
locations = doc.getElementsByTagName('location')
if locations.getLength > 1
fprintf('%d places called %s found\n', locations.getLength, place);
end
location = locations.item(0); % take the first return
assert(~isempty(location), 'location not found');
node_lat = location.getElementsByTagName('lat');
el = node_lat.item(0);
lat = str2num( el.getFirstChild.getData );
node_lon = location.getElementsByTagName('lng');
el = node_lon.item(0);
lon = str2num( el.getFirstChild.getData );
else
lat = args{1};
lon = args{2};
if length(args) == 3
zoom = args{3};
else
zoom = 18;
end
end
% now read the map
baseurl = sprintf('https://maps.googleapis.com/maps/api/staticmap?center=%.6f,%.6f&zoom=%d&size=%dx%d&scale=%d&format=png&maptype=%s&key=%s&sensor=false', lat, lon, zoom, ev.width, ev.height, opt.scale, opt.type, ev.key);
opturl = '';
if ~opt.roadnames
opturl = strcat(opturl, '&style=feature:road|element:labels|visibility:off');
end
if ~opt.placenames
opturl = strcat(opturl, '&style=feature:administrative|element:labels.text|visibility:off&style=feature:poi|visibility:off');
end
if opt.onlyroads
opturl = strcat(opturl, ['&style=feature:landscape|element:geometry.fill|color:0x000000|visibility:on'...
'&style=feature:landscape|element:labels|visibility:off'...
'&style=feature:administrative|visibility:off'...
'&style=feature:road|element:geometry.fill|color:0xffffff|visibility:on'...
'&style=feature:road|element:labels|visibility:off'...
'&style=feature:poi|element:all|visibility:off'...
'&style=feature:transit|element:all|visibility:off'...
'&style=feature:water|element:all|visibility:off'...
]);
end
%opturl = urlencode(opturl)
opturl = strrep(opturl, '|', '%7C');
%opturl
try
[idx,cmap] = imread([baseurl opturl], 'png');
catch
error('Failed to read the URL: check your Google API key');
end
cmap = iint(cmap);
% apply the color map
view = cmap(idx(:)+1,:);
% knock it into shape
view = shiftdim( reshape(view', [3 size(idx)]), 1);
view = ev.convert(view);
if nargout == 0
idisp(view);
else
im = view;
end
if nargout > 1
% compute the northing/easting at each pixel.
%
% the following lines of code from get_google_map by Val Schmidt
% ESTIMATE BOUNDS OF IMAGE:
%
% Normally one must specify a center (lat,lon) and zoom level for a map.
% Zoom Notes:
% ZL: 15, hxw = 640x640, image dim: 2224.91 x 2224.21 (mean 2224.56)
% ZL: 16, hxw = 640x640, image dim: 1128.01m x 1111.25m (mean 1119.63)
% This gives an equation of roughly (ZL-15)*3.4759 m/pix * pixels
% So for an image at ZL 16, the LHS bound is
% LHS = centerlonineastings - (zl-15) * 3.4759 * 640/2;
[lonutm latutm zone] = deg2utm(lat,lon);
Hdim = (2^15/2^zoom) * 3.4759 * ev.width;
Vdim = (2^15/2^zoom) * 3.4759 * ev.height;
ell = lonutm - Hdim/2;
nll = latutm - Vdim/2;
eur = lonutm + Hdim/2;
nur = latutm + Vdim/2;
Nvec = linspace(nur,nll,ev.height); % lat is highest at image row 1
Evec = linspace(ell,eur,ev.width);
[E,N] = meshgrid(Evec, Nvec);
end
end
function close()
end
function s = char(m)
%EarthView.char Convert to string
%
% EV.char() is a string representing the state of the EarthView object in
% human readable form.
%
% See also EarthView.display.
s = sprintf('EarthView: %d x %d', ...
m.width, m.height);
s2 = char@ImageSource(m);
if ~isempty(s2)
s = strvcat(s, strcat(' - ', s2));
end
end
function paramSet(varargin)
end
end
end