1
+ import os
1
2
import numpy as np
2
3
from shapefile import Reader
3
4
4
5
lsd = 5
5
6
7
+ UTILS_DIR = os .path .dirname (os .path .abspath (__file__ ))
8
+ OUTPUT_DIR = os .path .join (UTILS_DIR , '..' , 'lib' , 'mpl_toolkits' ,
9
+ 'basemap' , 'data' )
10
+
11
+ # Folder where GSHHG shapefiles were extracted. Change if needed
12
+ GSHHS_DIR = UTILS_DIR
13
+
6
14
def quantize (data ,least_significant_digit ):
7
15
"""
8
16
quantize data to improve compression. data is quantized using
@@ -25,8 +33,8 @@ def quantize(data,least_significant_digit):
25
33
def get_coast_polygons (resolution ):
26
34
polymeta = []; polybounds = []
27
35
for level in [1 ,2 ,3 ,5 ]:
28
- filename = 'GSHHS_shp/%s/GSHHS_%s_L%s' % ( resolution , resolution , level )
29
- #filename = 'WDBII_shp/%s/WDBII_border_%s_L%s' % (resolution, resolution, level)
36
+ filename = os . path . join ( GSHHS_DIR , 'GSHHS_shp/' , resolution ,
37
+ 'GSHHS_{}_L{}' . format (resolution , level ) )
30
38
print filename
31
39
shf = Reader (filename )
32
40
fields = shf .fields
@@ -59,9 +67,11 @@ def get_coast_polygons(resolution):
59
67
def get_wdb_boundaries (resolution ,level ,rivers = False ):
60
68
polymeta = []; polybounds = []
61
69
if rivers :
62
- filename = 'WDBII_shp/%s/WDBII_river_%s_L%02i' % (resolution , resolution , level )
70
+ filename = os .path .join (GSHHS_DIR , 'WDBII_shp' , resolution ,
71
+ 'WDBII_river_{}_L{:02}' .format (resolution , level ))
63
72
else :
64
- filename = 'WDBII_shp/%s/WDBII_border_%s_L%s' % (resolution , resolution , level )
73
+ filename = os .path .join (GSHHS_DIR , 'WDBII_shp' , resolution ,
74
+ 'WDBII_border_{}_L{}' .format (resolution , level ))
65
75
print filename
66
76
shf = Reader (filename )
67
77
fields = shf .fields
@@ -90,8 +100,8 @@ def get_wdb_boundaries(resolution,level,rivers=False):
90
100
# read in coastline data (only those polygons whose area > area_thresh).
91
101
for resolution in ['c' ,'l' ,'i' ,'h' ,'f' ]:
92
102
poly , polymeta = get_coast_polygons (resolution )
93
- f = open ('gshhs_' + resolution + '.dat' , 'wb' )
94
- f2 = open ('gshhsmeta_' + resolution + '.dat' , 'w' )
103
+ f = open (os . path . join ( OUTPUT_DIR , 'gshhs_' + resolution + '.dat' ), 'wb' )
104
+ f2 = open (os . path . join ( OUTPUT_DIR , 'gshhsmeta_' + resolution + '.dat' ), 'w' )
95
105
offset = 0
96
106
for p ,pm in zip (poly ,polymeta ):
97
107
typ = pm [0 ]; area = pm [1 ]; south = pm [2 ]; north = pm [3 ]; npts = pm [4 ]
@@ -106,8 +116,8 @@ def get_wdb_boundaries(resolution,level,rivers=False):
106
116
107
117
for resolution in ['c' ,'l' ,'i' ,'h' ,'f' ]:
108
118
poly , polymeta = get_wdb_boundaries (resolution ,1 )
109
- f = open ('countries_' + resolution + '.dat' , 'wb' )
110
- f2 = open ('countriesmeta_' + resolution + '.dat' , 'w' )
119
+ f = open (os . path . join ( OUTPUT_DIR , 'countries_' + resolution + '.dat' ), 'wb' )
120
+ f2 = open (os . path . join ( OUTPUT_DIR , 'countriesmeta_' + resolution + '.dat' ), 'w' )
111
121
offset = 0
112
122
for p ,pm in zip (poly ,polymeta ):
113
123
typ = pm [0 ]; area = pm [1 ]; south = pm [2 ]; north = pm [3 ]; npts = pm [4 ]
@@ -122,8 +132,8 @@ def get_wdb_boundaries(resolution,level,rivers=False):
122
132
123
133
for resolution in ['c' ,'l' ,'i' ,'h' ,'f' ]:
124
134
poly , polymeta = get_wdb_boundaries (resolution ,2 )
125
- f = open ('states_' + resolution + '.dat' , 'wb' )
126
- f2 = open ('statesmeta_' + resolution + '.dat' , 'w' )
135
+ f = open (os . path . join ( OUTPUT_DIR , 'states_' + resolution + '.dat' ), 'wb' )
136
+ f2 = open (os . path . join ( OUTPUT_DIR , 'statesmeta_' + resolution + '.dat' ), 'w' )
127
137
offset = 0
128
138
for p ,pm in zip (poly ,polymeta ):
129
139
typ = pm [0 ]; area = pm [1 ]; south = pm [2 ]; north = pm [3 ]; npts = pm [4 ]
@@ -137,8 +147,8 @@ def get_wdb_boundaries(resolution,level,rivers=False):
137
147
f2 .close ()
138
148
139
149
for resolution in ['c' ,'l' ,'i' ,'h' ,'f' ]:
140
- f = open ('rivers_' + resolution + '.dat' , 'wb' )
141
- f2 = open ('riversmeta_' + resolution + '.dat' , 'w' )
150
+ f = open (os . path . join ( OUTPUT_DIR , 'rivers_' + resolution + '.dat' ), 'wb' )
151
+ f2 = open (os . path . join ( OUTPUT_DIR , 'riversmeta_' + resolution + '.dat' ), 'w' )
142
152
for level in range (1 ,12 ):
143
153
poly , polymeta = get_wdb_boundaries (resolution ,level ,rivers = True )
144
154
offset = 0
0 commit comments