-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcartopy_mapping.py
163 lines (141 loc) · 6.51 KB
/
cartopy_mapping.py
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
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy
import matplotlib.pyplot as plt
import pyproj
import os, glob
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib import patheffects
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import numpy as np
from owslib.wmts import TileMatrixSetLink, TileMatrixLimits, _TILE_MATRIX_SET_TAG, _TILE_MATRIX_SET_LIMITS_TAG, _TILE_MATRIX_LIMITS_TAG
def custom_from_elements(link_elements):
links = []
for link_element in link_elements:
matrix_set_elements = link_element.findall(_TILE_MATRIX_SET_TAG)
if len(matrix_set_elements) == 0:
raise ValueError('Missing TileMatrixSet in %s' % link_element)
elif len(matrix_set_elements) > 1:
set_limits_elements = link_element.findall(
_TILE_MATRIX_SET_LIMITS_TAG)
if set_limits_elements:
raise ValueError('Multiple instances of TileMatrixSet'
' plus TileMatrixSetLimits in %s' %
link_element)
for matrix_set_element in matrix_set_elements:
uri = matrix_set_element.text.strip()
links.append(TileMatrixSetLink(uri))
else:
uri = matrix_set_elements[0].text.strip()
tilematrixlimits = {}
path = '%s/%s' % (_TILE_MATRIX_SET_LIMITS_TAG,
_TILE_MATRIX_LIMITS_TAG)
for limits_element in link_element.findall(path):
tml = TileMatrixLimits(limits_element)
if tml.tilematrix:
tilematrixlimits[tml.tilematrix] = tml
links.append(TileMatrixSetLink(uri, tilematrixlimits))
return links
TileMatrixSetLink.from_elements = custom_from_elements
def plot_map(df,show=True):
# make figure look nice
extent = (-139.2, -139.1, 61.325, 61.385)
fig = plt.figure(figsize=(15,15))
ax = fig.add_axes([0, 0, 1, 1],projection=ccrs.PlateCarree())
ax.set_extent(extent, ccrs.Geodetic())
# add color shaded relief
url = 'https://map1.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
layer = 'ASTER_GDEM_Color_Shaded_Relief'
ax.add_wmts(url, layer,alpha=0.8)
# plot data
plt.scatter(df['Long'],df['Lat'],marker='o',
s=300,color='blue',edgecolor='black',label='Station',zorder=4)
# plot station names
x_label_offset=0.00125
y_label_offset=-0.0025
for x, y, label in zip(df['Long'],df['Lat'], df['Station']):
ax.text(x+x_label_offset, y+y_label_offset, label,
ha='left',va='bottom',color='white', size=15)
# add lat-lon lines.
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=0.8, color='black')
gl.xlocator = mticker.FixedLocator(np.linspace(extent[0],extent[1],num=5))
gl.ylocator = mticker.FixedLocator(np.linspace(extent[2],extent[3],num=5))
gl.xlabels_bottom=True
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylabels_left=True
gl.xlabel_style = {'color': 'black','size': 15}
gl.ylabel_style = {'color': 'black','size': 15}
if show:
plt.show()
return ax
def plot_map_mt_objs(mt_objs,show=True):
# parse data objects
coords = [(mt_obj.lat, mt_obj.lon ,mt_obj.station) for mt_obj in mt_objs]
coords = np.asarray(coords,dtype=np.float64).T
extent = (-139.2, -139.1, 61.325, 61.385)
fig = plt.figure(figsize=(15,15))
ax = fig.add_axes([0, 0, 1, 1],projection=ccrs.PlateCarree())
ax.set_extent(extent, ccrs.Geodetic())
# add color shaded relief
url = 'https://map1.vis.earthdata.nasa.gov/wmts-geo/wmts.cgi'
layer = 'ASTER_GDEM_Color_Shaded_Relief'
ax.add_wmts(url, layer,alpha=0.8)
# plot data
plt.scatter(coords[1,:],coords[0,:],marker='o',
s=300,color='blue',edgecolor='black',label='Station',zorder=4)
# plot station names
x_label_offset=0.00125
y_label_offset=-0.0025
for x, y, label in zip(coords[1,:],coords[0,:],coords[2,:]):
str_label = str(int(float(label)))
x_coord = float(x) + x_label_offset
y_coord = float(y) + y_label_offset
ax.text(x_coord, y_coord, str_label,ha='left',va='bottom',color='white', size=15)
# add lat-lon lines.
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=0.8, color='black')
gl.xlocator = mticker.FixedLocator(np.linspace(extent[0],extent[1],num=5))
gl.ylocator = mticker.FixedLocator(np.linspace(extent[2],extent[3],num=5))
gl.xlabels_bottom=True
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.ylabels_left=True
gl.xlabel_style = {'color': 'black','size': 15}
gl.ylabel_style = {'color': 'black','size': 15}
if show:
plt.show()
return ax
def convert_wgs_to_utm(origin):
# from https://stackoverflow.com/questions/40132542/get-a-cartesian-projection-accurate-around-a-lat-lng-pair
utm_band = str(int((np.floor((origin[0] + 180) / 6 ) % 60) + 1))
if len(utm_band) == 1:
utm_band = '0'+utm_band
if origin[1] >= 0:
epsg_code = '326' + utm_band
else:
epsg_code = '327' + utm_band
return epsg_code
def plot_induction_arrows(ax,df,color='white',scale=1000,show=True,zorder=4):
this_df = df.copy()
this_df = this_df.dropna(axis=0).reset_index()
SCALE = scale
origin = (this_df['Long'].values[0], this_df['Lat'].values[0])
utm_code = convert_wgs_to_utm(origin)
crs_wgs = pyproj.Proj(init='epsg:4326',) # assuming you're using WGS84 geographic
crs_utm = pyproj.Proj(init='epsg:{0}'.format(utm_code))
east_0, north_0 = pyproj.transform(crs_wgs, crs_utm,
this_df['Long'].values, this_df['Lat'].values,
always_xy='True')
east, north = pyproj.transform(crs_utm, crs_wgs,
this_df['East'].values*SCALE +east_0, this_df['North'].values*SCALE +north_0,
always_xy='True')
east = np.asarray(east) - this_df['Long']
north= np.asarray(north) - this_df['Lat']
for index, row in this_df.iterrows():
ax.arrow(row['Long'],row['Lat'],east[index],north[index],zorder=zorder,
facecolor=color,head_width=0.001,head_length=0.002,width=0.0006)
if show:
plt.show()