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
|
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) 2012-2019 Satpy developers
#
# This file is part of satpy.
#
# satpy is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your option) any later
# version.
#
# satpy is distributed in the hope that it will be useful, but WITHOUT ANY
# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
# A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with
# satpy. If not, see <http://www.gnu.org/licenses/>.
"""Convert human coordinates (lon and lat) to an area definition.
Here is a usage example.
python coord2area_def.py france stere 42.0 51.5 -5.5 8.0 1.5
The arguments are "name proj min_lat max_lat min_lon max_lon resolution(km)".
The command above yelds the following result.
### +proj=stere +lat_0=46.75 +lon_0=1.25 +ellps=WGS84
france:
description: france
projection:
proj: stere
ellps: WGS84
lat_0: 46.75
lon_0: 1.25
shape:
height: 703
width: 746
area_extent:
lower_left_xy: [-559750.381098, -505020.675776]
upper_right_xy: [559750.381098, 549517.351948]
The first commented line is just a sum-up. The value of "description" can be changed to any descriptive text.
Such a custom yaml configuration can be profitably saved in a local areas.yaml configuration file that won't be
overridden by future updates of SatPy package. For that purpose the local processing script may have suitable
lines as reported below.
# set PPP_CONFIG_DIR for custom composites
import os
os.environ['PPP_CONFIG_DIR'] = '/my_local_path/for_satpy_configuration'
As a further functionality this script may give a quick display of the defined area,
provided the path for the GSHHG library is supplied via the "-s" option
and the modules PyCoast, Pillow and AggDraw have been installed.
python coord2area_def.py france stere 42.0 51.5 -5.5 8.0 1.5 -s /path/for/gshhs/library
The command above would first print the seen area definition and then launch a casual representation
of the area relying on the information about borders involved.
"""
import argparse
import sys
from pyproj import Proj
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("name",
help="The name of the area.")
parser.add_argument("proj",
help="The projection to use. Use proj.4 names, like 'stere', 'merc'...")
parser.add_argument("min_lat",
help="The the latitude of the bottom of the area",
type=float)
parser.add_argument("max_lat",
help="The the latitude of the top of the area",
type=float)
parser.add_argument("min_lon",
help="The the longitude of the left of the area",
type=float)
parser.add_argument("max_lon",
help="The the longitude of the right of the area",
type=float)
parser.add_argument("resolution",
help="The resolution of the area (in km)",
type=float)
parser.add_argument("-s", "--shapes",
help="Show a preview of the area using the coastlines in this directory")
args = parser.parse_args()
name = args.name
proj = args.proj
left = args.min_lon
right = args.max_lon
up = args.min_lat
down = args.max_lat
res = args.resolution * 1000
lat_0 = (up + down) / 2
lon_0 = (right + left) / 2
p = Proj(proj=proj, lat_0=lat_0, lon_0=lon_0, ellps="WGS84")
left_ex1, up_ex1 = p(left, up)
right_ex1, up_ex2 = p(right, up)
left_ex2, down_ex1 = p(left, down)
right_ex2, down_ex2 = p(right, down)
left_ex3, dummy = p(left, lat_0)
right_ex3, dummy = p(right, lat_0)
area_extent = (min(left_ex1, left_ex2, left_ex3),
min(up_ex1, up_ex2),
max(right_ex1, right_ex2, right_ex3),
max(down_ex1, down_ex2))
xsize = int(round((area_extent[2] - area_extent[0]) / res))
ysize = int(round((area_extent[3] - area_extent[1]) / res))
proj4_string = "+" + \
" +".join(("proj=" + proj + ",lat_0=" + str(lat_0) +
",lon_0=" + str(lon_0) + ",ellps=WGS84").split(","))
print("### " + proj4_string)
print()
print(name + ":")
print(" description: " + name)
print(" projection:")
print(" proj: " + proj)
print(" ellps: WGS84")
print(" lat_0: " + str(lat_0))
print(" lon_0: " + str(lon_0))
print(" shape:")
print(" height: " + str(ysize))
print(" width: " + str(xsize))
print(" area_extent:")
print(" lower_left_xy: [%f, %f]" % (area_extent[0], area_extent[1]))
print(" upper_right_xy: [%f, %f]" % (area_extent[2], area_extent[3]))
if args.shapes is None:
sys.exit(0)
from PIL import Image
from pycoast import ContourWriterAGG
img = Image.new("RGB", (xsize, ysize))
area_def = (proj4_string, area_extent)
cw = ContourWriterAGG(args.shapes)
cw.add_coastlines(img, (proj4_string, area_extent),
resolution="l", width=0.5)
cw.add_grid(img, area_def, (10.0, 10.0), (2.0, 2.0), write_text=False, outline="white", outline_opacity=175,
width=1.0, minor_outline="white", minor_outline_opacity=175, minor_width=0.2, minor_is_tick=False)
img.show()
|
