New upstream version 2.1.3+ds
Bas Couwenberg
5 years ago
79 | 79 | |
80 | 80 | script: |
81 | 81 | - python -c "import pyproj; pyproj.Proj(init='epsg:4269')" |
82 | - nose2 -v | |
83 | - nose2 --coverage pyproj --coverage-report xml # coverage report on separate line as does not show failures | |
82 | - py.test --cov-report term-missing --cov=pyproj -v -s | |
84 | 83 | |
85 | 84 | after_success: |
86 | 85 | - coveralls |
116 | 116 | test_script: |
117 | 117 | # Run the project tests |
118 | 118 | - "%CMD_IN_ENV% python -c \"import pyproj; pyproj.Proj(init='epsg:4269')\"" |
119 | - "%CMD_IN_ENV% python unittest/test.py -v" | |
119 | - "%CMD_IN_ENV% py.test --cov-report term-missing --cov=pyproj -v -s" | |
120 | 120 | |
121 | 121 | after_test: |
122 | 122 | # If tests are successful, create binary packages for the project. |
0 | [unittest] | |
1 | start-dir = unittest | |
2 | ||
3 | [coverage] | |
4 | always-on = True | |
5 | coverage = pyproj/_proj.pyx | |
6 | pyproj | |
7 | coverage-report = term |
46 | 46 | LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, |
47 | 47 | NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN |
48 | 48 | CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. """ |
49 | __version__ = "2.1.2" | |
49 | __version__ = "2.1.3" | |
50 | 50 | __all__ = [ |
51 | 51 | "Proj", |
52 | 52 | "Geod", |
78 | 78 | pyproj.transformer, verbose=verbose |
79 | 79 | ) |
80 | 80 | |
81 | return failure_count + failure_count_crs + failure_count_geod | |
81 | return failure_count + failure_count_crs + failure_count_geod + failure_count_transform | |
82 | 82 | |
83 | 83 | |
84 | 84 | if __name__ == "__main__": |
350 | 350 | """ |
351 | 351 | Returns |
352 | 352 | ------- |
353 | Ellipsoid: The CRS ellipsoid object with associated attributes. | |
353 | Ellipsoid: The ellipsoid object with associated attributes. | |
354 | 354 | """ |
355 | 355 | if self._ellipsoid is not None: |
356 | 356 | return self._ellipsoid |
362 | 362 | """ |
363 | 363 | Returns |
364 | 364 | ------- |
365 | AreaOfUse: The CRS area of use object with associated attributes. | |
365 | AreaOfUse: The area of use object with associated attributes. | |
366 | 366 | """ |
367 | 367 | if self._area_of_use is not None: |
368 | 368 | return self._area_of_use |
386 | 386 | """ |
387 | 387 | Returns |
388 | 388 | ------- |
389 | :obj:`pyproj.CRS`: The datum as a CRS. | |
389 | pyproj.CRS: The datum as a CRS. | |
390 | 390 | """ |
391 | 391 | cdef PJ *projobj = proj_crs_get_datum(self.projctx, self.projobj) |
392 | 392 | if projobj == NULL: |
435 | 435 | return None |
436 | 436 | |
437 | 437 | def to_geodetic(self): |
438 | """return a new CRS instance which is the geographic (lat/lon) | |
439 | coordinate version of the current projection""" | |
438 | """ | |
439 | ||
440 | Returns | |
441 | ------- | |
442 | pyproj.CRS: The geographic (lat/lon) CRS from the current CRS. | |
443 | ||
444 | """ | |
440 | 445 | cdef PJ * projobj |
441 | 446 | projobj = proj_crs_get_geodetic_crs(self.projctx, self.projobj) |
442 | 447 | if projobj == NULL: |
449 | 454 | def to_wkt(self, version="WKT2_2018"): |
450 | 455 | """ |
451 | 456 | Convert the projection to a WKT string. |
457 | ||
458 | Version options: | |
459 | - WKT2_2015 | |
460 | - WKT2_2015_SIMPLIFIED | |
461 | - WKT2_2018 | |
462 | - WKT2_2018_SIMPLIFIED | |
463 | - WKT1_GDAL | |
464 | - WKT1_ESRI | |
465 | ||
452 | 466 | |
453 | 467 | Parameters |
454 | 468 | ---------- |
455 | 469 | version: str |
456 | 470 | The version of the WKT output. Default is WKT2_2018. |
457 | ||
471 | ||
458 | 472 | Returns |
459 | 473 | ------- |
460 | 474 | str: The WKT string. |
0 | import os | |
0 | 1 | |
1 | 2 | from libc.stdlib cimport malloc, free |
2 | 3 | |
14 | 15 | |
15 | 16 | cdef PJ_CONTEXT* get_pyproj_context(): |
16 | 17 | data_dir = get_data_dir() |
17 | data_dir_list = data_dir.split(";") | |
18 | data_dir_list = data_dir.split(os.pathsep) | |
18 | 19 | cdef PJ_CONTEXT* pyproj_context = NULL |
19 | 20 | cdef char **c_data_dir = <char **>malloc(len(data_dir_list) * sizeof(char*)) |
20 | 21 | try: |
90 | 90 | if errcheck: |
91 | 91 | err = proj_errno(self.projpj) |
92 | 92 | if err != 0: |
93 | raise ProjError(proj_errno_string(err)) | |
93 | raise ProjError(pystrdecode(proj_errno_string(err))) | |
94 | 94 | # since HUGE_VAL can be 'inf', |
95 | 95 | # change it to a real (but very large) number. |
96 | 96 | # also check for NaNs. |
148 | 148 | if errcheck: |
149 | 149 | err = proj_errno(self.projpj) |
150 | 150 | if err != 0: |
151 | raise ProjError(proj_errno_string(err)) | |
151 | raise ProjError(pystrdecode(proj_errno_string(err))) | |
152 | 152 | # since HUGE_VAL can be 'inf', |
153 | 153 | # change it to a real (but very large) number. |
154 | 154 | # also check for NaNs. |
38 | 38 | cdef _Transformer transformer = _Transformer() |
39 | 39 | transformer.projpj = proj_create_crs_to_crs( |
40 | 40 | transformer.projctx, |
41 | _Transformer._definition_from_object(proj_from), | |
42 | _Transformer._definition_from_object(proj_to), | |
41 | cstrencode(proj_from.crs.srs), | |
42 | cstrencode(proj_to.crs.srs), | |
43 | 43 | NULL) |
44 | 44 | if transformer.projpj is NULL: |
45 | 45 | raise ProjError("Error creating CRS to CRS.") |
46 | 46 | transformer.set_radians_io() |
47 | transformer.projections_exact_same = proj_from.is_exact_same(proj_to) | |
48 | transformer.projections_equivalent = proj_from == proj_to | |
47 | transformer.projections_exact_same = proj_from.crs.is_exact_same(proj_to.crs) | |
48 | transformer.projections_equivalent = proj_from.crs == proj_to.crs | |
49 | 49 | transformer.skip_equivalent = skip_equivalent |
50 | 50 | transformer.is_pipeline = False |
51 | 51 | return transformer |
62 | 62 | return transformer |
63 | 63 | |
64 | 64 | @staticmethod |
65 | def from_crs(crs_from, crs_to, skip_equivalent=False): | |
66 | if not isinstance(crs_from, CRS): | |
67 | crs_from = CRS.from_user_input(crs_from) | |
68 | if not isinstance(crs_to, CRS): | |
69 | crs_to = CRS.from_user_input(crs_to) | |
70 | transformer = _Transformer._init_crs_to_crs(crs_from, crs_to, skip_equivalent=skip_equivalent) | |
71 | transformer.input_geographic = crs_from.is_geographic | |
72 | transformer.output_geographic = crs_to.is_geographic | |
73 | return transformer | |
74 | ||
75 | @staticmethod | |
76 | 65 | def from_pipeline(const char *proj_pipeline): |
77 | 66 | cdef _Transformer transformer = _Transformer() |
78 | 67 | |
84 | 73 | transformer.is_pipeline = True |
85 | 74 | return transformer |
86 | 75 | |
87 | @staticmethod | |
88 | def _definition_from_object(in_proj): | |
89 | """ | |
90 | Parameters | |
91 | ---------- | |
92 | in_proj: :obj:`pyproj.Proj` or :obj:`pyproj.CRS` | |
93 | ||
94 | Returns | |
95 | ------- | |
96 | char*: Definition string for `proj_create_crs_to_crs`. | |
97 | ||
98 | """ | |
99 | if isinstance(in_proj, Proj): | |
100 | return cstrencode(in_proj.crs.srs) | |
101 | return cstrencode(in_proj.srs) | |
102 | ||
103 | def _transform(self, inx, iny, inz, radians, errcheck=False): | |
76 | def _transform(self, inx, iny, inz, intime, radians, errcheck=False): | |
104 | 77 | if self.projections_exact_same or (self.projections_equivalent and self.skip_equivalent): |
105 | 78 | return |
106 | 79 | # private function to call pj_transform |
107 | 80 | cdef void *xdata |
108 | 81 | cdef void *ydata |
109 | 82 | cdef void *zdata |
83 | cdef void *tdata | |
110 | 84 | cdef double *xx |
111 | 85 | cdef double *yy |
112 | 86 | cdef double *zz |
113 | cdef Py_ssize_t buflenx, bufleny, buflenz, npts, i | |
87 | cdef double *tt | |
88 | cdef Py_ssize_t buflenx, bufleny, buflenz, buflent, npts, iii | |
114 | 89 | cdef int err |
115 | 90 | if PyObject_AsWriteBuffer(inx, &xdata, &buflenx) <> 0: |
116 | 91 | raise ProjError |
121 | 96 | raise ProjError |
122 | 97 | else: |
123 | 98 | buflenz = bufleny |
124 | if not (buflenx == bufleny == buflenz): | |
125 | raise ProjError('x,y and z must be same size') | |
99 | if intime is not None: | |
100 | if PyObject_AsWriteBuffer(intime, &tdata, &buflent) <> 0: | |
101 | raise ProjError | |
102 | else: | |
103 | buflent = bufleny | |
104 | ||
105 | if not buflenx or not (buflenx == bufleny == buflenz == buflent): | |
106 | raise ProjError('x,y,z, and time must be same size') | |
126 | 107 | xx = <double *>xdata |
127 | 108 | yy = <double *>ydata |
128 | 109 | if inz is not None: |
129 | 110 | zz = <double *>zdata |
130 | 111 | else: |
131 | 112 | zz = NULL |
113 | if intime is not None: | |
114 | tt = <double *>tdata | |
115 | else: | |
116 | tt = NULL | |
132 | 117 | npts = buflenx//8 |
133 | 118 | |
134 | 119 | # degrees to radians |
135 | 120 | if not self.is_pipeline and not radians and self.input_radians: |
136 | for i from 0 <= i < npts: | |
137 | xx[i] = xx[i]*_DG2RAD | |
138 | yy[i] = yy[i]*_DG2RAD | |
121 | for iii from 0 <= iii < npts: | |
122 | xx[iii] = xx[iii]*_DG2RAD | |
123 | yy[iii] = yy[iii]*_DG2RAD | |
139 | 124 | # radians to degrees |
140 | 125 | elif not self.is_pipeline and radians and not self.input_radians and self.input_geographic: |
141 | for i from 0 <= i < npts: | |
142 | xx[i] = xx[i]*_RAD2DG | |
143 | yy[i] = yy[i]*_RAD2DG | |
144 | ||
145 | proj_trans_generic( | |
126 | for iii from 0 <= iii < npts: | |
127 | xx[iii] = xx[iii]*_RAD2DG | |
128 | yy[iii] = yy[iii]*_RAD2DG | |
129 | ||
130 | cdef int err_count = proj_trans_generic( | |
146 | 131 | self.projpj, |
147 | 132 | PJ_FWD, |
148 | 133 | xx, _DOUBLESIZE, npts, |
149 | 134 | yy, _DOUBLESIZE, npts, |
150 | 135 | zz, _DOUBLESIZE, npts, |
151 | NULL, 0, 0, | |
136 | tt, _DOUBLESIZE, npts, | |
152 | 137 | ) |
153 | 138 | cdef int errno = proj_errno(self.projpj) |
154 | 139 | if errno and errcheck: |
155 | 140 | raise ProjError("proj_trans_generic error: {}".format( |
156 | 141 | pystrdecode(proj_errno_string(errno)))) |
142 | elif err_count and errcheck: | |
143 | raise ProjError("{} proj_trans_generic error(s)".format(err_count)) | |
157 | 144 | |
158 | 145 | # radians to degrees |
159 | 146 | if not self.is_pipeline and not radians and self.output_radians: |
160 | for i from 0 <= i < npts: | |
161 | xx[i] = xx[i]*_RAD2DG | |
162 | yy[i] = yy[i]*_RAD2DG | |
147 | for iii from 0 <= iii < npts: | |
148 | xx[iii] = xx[iii]*_RAD2DG | |
149 | yy[iii] = yy[iii]*_RAD2DG | |
163 | 150 | # degrees to radians |
164 | 151 | elif not self.is_pipeline and radians and not self.output_radians and self.output_geographic: |
165 | for i from 0 <= i < npts: | |
166 | xx[i] = xx[i]*_DG2RAD | |
167 | yy[i] = yy[i]*_DG2RAD | |
152 | for iii from 0 <= iii < npts: | |
153 | xx[iii] = xx[iii]*_DG2RAD | |
154 | yy[iii] = yy[iii]*_DG2RAD | |
168 | 155 | |
169 | 156 | |
170 | 157 | def _transform_sequence(self, Py_ssize_t stride, inseq, bint switch, |
171 | radians, errcheck=False): | |
158 | time_3rd, radians, errcheck=False): | |
172 | 159 | if self.projections_exact_same or (self.projections_equivalent and self.skip_equivalent): |
173 | 160 | return |
174 | 161 | # private function to itransform function |
178 | 165 | double *x |
179 | 166 | double *y |
180 | 167 | double *z |
181 | Py_ssize_t buflen, npts, i, j | |
168 | double *tt | |
169 | Py_ssize_t buflen, npts, iii, jjj | |
182 | 170 | int err |
183 | 171 | |
184 | 172 | if stride < 2: |
191 | 179 | |
192 | 180 | # degrees to radians |
193 | 181 | if not self.is_pipeline and not radians and self.input_radians: |
194 | for i from 0 <= i < npts: | |
195 | j = stride*i | |
196 | coords[j] *= _DG2RAD | |
197 | coords[j+1] *= _DG2RAD | |
182 | for iii from 0 <= iii < npts: | |
183 | jjj = stride*iii | |
184 | coords[jjj] *= _DG2RAD | |
185 | coords[jjj+1] *= _DG2RAD | |
198 | 186 | # radians to degrees |
199 | 187 | elif not self.is_pipeline and radians and not self.input_radians and self.input_geographic: |
200 | for i from 0 <= i < npts: | |
201 | j = stride*i | |
202 | coords[j] *= _RAD2DG | |
203 | coords[j+1] *= _RAD2DG | |
188 | for iii from 0 <= iii < npts: | |
189 | jjj = stride*iii | |
190 | coords[jjj] *= _RAD2DG | |
191 | coords[jjj+1] *= _RAD2DG | |
204 | 192 | |
205 | 193 | if not switch: |
206 | 194 | x = coords |
209 | 197 | x = coords + 1 |
210 | 198 | y = coords |
211 | 199 | |
212 | if stride == 2: | |
200 | # z coordinate | |
201 | if stride == 4 or (stride == 3 and not time_3rd): | |
202 | z = coords + 2 | |
203 | else: | |
213 | 204 | z = NULL |
214 | else: | |
215 | z = coords + 2 | |
216 | ||
217 | proj_trans_generic ( | |
205 | # time | |
206 | if stride == 3 and time_3rd: | |
207 | tt = coords + 2 | |
208 | elif stride == 4: | |
209 | tt = coords + 3 | |
210 | else: | |
211 | tt = NULL | |
212 | ||
213 | cdef int err_count = proj_trans_generic ( | |
218 | 214 | self.projpj, |
219 | 215 | PJ_FWD, |
220 | 216 | x, stride*_DOUBLESIZE, npts, |
221 | 217 | y, stride*_DOUBLESIZE, npts, |
222 | 218 | z, stride*_DOUBLESIZE, npts, |
223 | NULL, 0, 0, | |
219 | tt, stride*_DOUBLESIZE, npts, | |
224 | 220 | ) |
225 | ||
226 | 221 | cdef int errno = proj_errno(self.projpj) |
227 | 222 | if errno and errcheck: |
228 | 223 | raise ProjError("proj_trans_generic error: {}".format( |
229 | proj_errno_string(errno))) | |
224 | pystrdecode(proj_errno_string(errno)))) | |
225 | elif err_count and errcheck: | |
226 | raise ProjError("{} proj_trans_generic error(s)".format(err_count)) | |
227 | ||
230 | 228 | |
231 | 229 | # radians to degrees |
232 | 230 | if not self.is_pipeline and not radians and self.output_radians: |
233 | for i from 0 <= i < npts: | |
234 | j = stride*i | |
235 | coords[j] *= _RAD2DG | |
236 | coords[j+1] *= _RAD2DG | |
231 | for iii from 0 <= iii < npts: | |
232 | jjj = stride*iii | |
233 | coords[jjj] *= _RAD2DG | |
234 | coords[jjj+1] *= _RAD2DG | |
237 | 235 | # degrees to radians |
238 | 236 | elif not self.is_pipeline and radians and not self.output_radians and self.output_geographic: |
239 | for i from 0 <= i < npts: | |
240 | j = stride*i | |
241 | coords[j] *= _DG2RAD | |
242 | coords[j+1] *= _DG2RAD | |
237 | for iii from 0 <= iii < npts: | |
238 | jjj = stride*iii | |
239 | coords[jjj] *= _DG2RAD | |
240 | coords[jjj+1] *= _DG2RAD |
159 | 159 | |
160 | 160 | Returns |
161 | 161 | ------- |
162 | CRS | |
162 | ~CRS | |
163 | 163 | """ |
164 | 164 | if int(code) <= 0: |
165 | 165 | raise CRSError("EPSG codes are positive integers") |
176 | 176 | |
177 | 177 | Returns |
178 | 178 | ------- |
179 | CRS | |
179 | ~CRS | |
180 | 180 | """ |
181 | 181 | if not proj_string: |
182 | 182 | raise CRSError("CRS is empty or invalid: {!r}".format(proj_string)) |
213 | 213 | |
214 | 214 | Returns |
215 | 215 | ------- |
216 | CRS | |
216 | ~CRS | |
217 | 217 | """ |
218 | 218 | if isinstance(value, _CRS): |
219 | 219 | return value |
232 | 232 | """ |
233 | 233 | Returns |
234 | 234 | ------- |
235 | :obj:`pyproj.Geod`: Geod object based on the CRS.ellipsoid. | |
235 | ~pyproj.geod.Geod: Geod object based on the ellipsoid. | |
236 | 236 | """ |
237 | 237 | if self.ellipsoid is None or not self.ellipsoid.ellipsoid_loaded: |
238 | 238 | return None |
34 | 34 | proj_data_dir: str |
35 | 35 | The path to rhe PROJ.4 data directory. |
36 | 36 | """ |
37 | set_data_dir(";".join([get_data_dir(), proj_data_dir])) | |
37 | set_data_dir(os.pathsep.join([get_data_dir(), proj_data_dir])) | |
38 | 38 | |
39 | 39 | |
40 | 40 | def get_data_dir(): |
72 | 72 | def valid_data_dirs(potential_data_dirs): |
73 | 73 | if potential_data_dirs is None: |
74 | 74 | return False |
75 | for proj_data_dir in potential_data_dirs.split(";"): | |
75 | for proj_data_dir in potential_data_dirs.split(os.pathsep): | |
76 | 76 | if valid_data_dir(proj_data_dir): |
77 | 77 | return True |
78 | 78 | break |
350 | 350 | ------- |
351 | 351 | bool: True if projection in geographic (lon/lat) coordinates. |
352 | 352 | """ |
353 | warnings.warn("'is_latlong()' is deprecated. Please use 'crs.is_geographic'.") | |
353 | warnings.warn( | |
354 | "'is_latlong()' is deprecated and will be removed in version 2.2.0. " | |
355 | "Please use 'crs.is_geographic'." | |
356 | ) | |
354 | 357 | return self.crs.is_geographic |
355 | 358 | |
356 | 359 | def is_geocent(self): |
359 | 362 | ------- |
360 | 363 | bool: True if projection in geocentric (x/y) coordinates |
361 | 364 | """ |
362 | warnings.warn("'is_geocent()' is deprecated. Please use 'crs.is_geocent'.") | |
365 | warnings.warn( | |
366 | "'is_geocent()' is deprecated and will be removed in version 2.2.0. " | |
367 | "Please use 'crs.is_geocent'." | |
368 | ) | |
363 | 369 | return self.is_geocent |
364 | 370 | |
365 | 371 | def definition_string(self): |
16 | 16 | NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN |
17 | 17 | CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.""" |
18 | 18 | |
19 | import warnings | |
19 | 20 | from array import array |
20 | 21 | from itertools import chain, islice |
21 | 22 | |
40 | 41 | |
41 | 42 | @staticmethod |
42 | 43 | def from_proj(proj_from, proj_to, skip_equivalent=False): |
43 | """Make a Transformer from a :obj:`pyproj.Proj` or input used to create one. | |
44 | """Make a Transformer from a :obj:`~pyproj.proj.Proj` or input used to create one. | |
44 | 45 | |
45 | 46 | Parameters |
46 | 47 | ---------- |
47 | proj_from: :obj:`pyproj.Proj` or input used to create one | |
48 | proj_from: :obj:`~pyproj.proj.Proj` or input used to create one | |
48 | 49 | Projection of input data. |
49 | proj_from: :obj:`pyproj.Proj` or input used to create one | |
50 | proj_to: :obj:`~pyproj.proj.Proj` or input used to create one | |
50 | 51 | Projection of output data. |
51 | 52 | skip_equivalent: bool, optional |
52 | 53 | If true, will skip the transformation operation if input and output |
54 | 55 | |
55 | 56 | Returns |
56 | 57 | ------- |
57 | :obj:`pyproj.Transformer` | |
58 | :obj:`~Transformer` | |
58 | 59 | |
59 | 60 | """ |
60 | 61 | |
66 | 67 | |
67 | 68 | @staticmethod |
68 | 69 | def from_crs(crs_from, crs_to, skip_equivalent=False): |
69 | """Make a Transformer from a :obj:`pyproj.CRS` or input used to create one. | |
70 | """Make a Transformer from a :obj:`~pyproj.crs.CRS` or input used to create one. | |
71 | ||
72 | .. warning:: `from_crs` is deprecated and will be removed in 2.2.0. | |
73 | Please use :meth:`~Transformer.from_proj` instead. | |
70 | 74 | |
71 | 75 | Parameters |
72 | 76 | ---------- |
73 | proj_from: :obj:`pyproj.CRS` or input used to create one | |
77 | crs_from: ~pyproj.crs.CRS or input used to create one | |
74 | 78 | Projection of input data. |
75 | proj_from: :obj:`pyproj.CRS` or input used to create one | |
79 | crs_to: ~pyproj.crs.CRS or input used to create one | |
76 | 80 | Projection of output data. |
77 | 81 | skip_equivalent: bool, optional |
78 | 82 | If true, will skip the transformation operation if input and output |
80 | 84 | |
81 | 85 | Returns |
82 | 86 | ------- |
83 | :obj:`pyproj.Transformer` | |
84 | ||
85 | """ | |
86 | transformer = Transformer() | |
87 | transformer._transformer = _Transformer.from_crs( | |
88 | crs_from, crs_to, skip_equivalent | |
87 | :obj:`~Transformer` | |
88 | ||
89 | """ | |
90 | warnings.warn( | |
91 | "`from_crs` is deprecated and will be removed in 2.2.0. " | |
92 | "Please use `from_proj` instead." | |
89 | 93 | ) |
90 | return transformer | |
94 | return Transformer.from_proj(crs_from, crs_to, skip_equivalent=skip_equivalent) | |
91 | 95 | |
92 | 96 | @staticmethod |
93 | 97 | def from_pipeline(proj_pipeline): |
102 | 106 | |
103 | 107 | Returns |
104 | 108 | ------- |
105 | :obj:`pyproj.Transformer` | |
109 | ~Transformer | |
106 | 110 | |
107 | 111 | """ |
108 | 112 | transformer = Transformer() |
109 | 113 | transformer._transformer = _Transformer.from_pipeline(cstrencode(proj_pipeline)) |
110 | 114 | return transformer |
111 | 115 | |
112 | def transform(self, xx, yy, zz=None, radians=False, errcheck=False): | |
116 | def transform(self, xx, yy, zz=None, tt=None, radians=False, errcheck=False): | |
113 | 117 | """ |
114 | 118 | Transform points between two coordinate systems. |
115 | 119 | |
121 | 125 | Input y coordinate(s). |
122 | 126 | zz: scalar or array (numpy or python), optional |
123 | 127 | Input z coordinate(s). |
128 | tt: scalar or array (numpy or python), optional | |
129 | Input time coordinate(s). | |
124 | 130 | radians: boolean, optional |
125 | 131 | If True, will expect input data to be in radians and will return radians |
126 | 132 | if the projection is geographic. Default is False (degrees). Ignored for |
151 | 157 | >>> xpjr, ypjr, zpjr = transprojr.transform(xpj, ypj, zpj, radians=True) |
152 | 158 | >>> "%.3f %.3f %.3f" % (xpjr, ypjr, zpjr) |
153 | 159 | '-2704026.010 -4253051.810 3895878.820' |
154 | >>> transformer = Transformer.from_crs("epsg:4326", 4326, skip_equivalent=True) | |
160 | >>> transformer = Transformer.from_proj("epsg:4326", 4326, skip_equivalent=True) | |
155 | 161 | >>> xeq, yeq = transformer.transform(33, 98) |
156 | 162 | >>> "%.0f %.0f" % (xeq, yeq) |
157 | 163 | '33 98' |
164 | 170 | inz, zisfloat, zislist, zistuple = _copytobuffer(zz) |
165 | 171 | else: |
166 | 172 | inz = None |
173 | if tt is not None: | |
174 | intime, tisfloat, tislist, tistuple = _copytobuffer(tt) | |
175 | else: | |
176 | intime = None | |
167 | 177 | # call pj_transform. inx,iny,inz buffers modified in place. |
168 | self._transformer._transform(inx, iny, inz, radians, errcheck=errcheck) | |
178 | self._transformer._transform( | |
179 | inx, iny, inz=inz, intime=intime, radians=radians, errcheck=errcheck | |
180 | ) | |
169 | 181 | # if inputs were lists, tuples or floats, convert back. |
170 | 182 | outx = _convertback(xisfloat, xislist, xistuple, inx) |
171 | 183 | outy = _convertback(yisfloat, yislist, xistuple, iny) |
184 | return_data = (outx, outy) | |
172 | 185 | if inz is not None: |
173 | outz = _convertback(zisfloat, zislist, zistuple, inz) | |
174 | return outx, outy, outz | |
175 | else: | |
176 | return outx, outy | |
177 | ||
178 | def itransform(self, points, switch=False, radians=False, errcheck=False): | |
186 | return_data += (_convertback(zisfloat, zislist, zistuple, inz),) | |
187 | if intime is not None: | |
188 | return_data += (_convertback(tisfloat, tislist, tistuple, intime),) | |
189 | return return_data | |
190 | ||
191 | def itransform( | |
192 | self, points, switch=False, time_3rd=False, radians=False, errcheck=False | |
193 | ): | |
179 | 194 | """ |
180 | 195 | Iterator/generator version of the function pyproj.Transformer.transform. |
181 | 196 | |
187 | 202 | switch: boolean, optional |
188 | 203 | If True x, y or lon,lat coordinates of points are switched to y, x |
189 | 204 | or lat, lon. Default is False. |
205 | time_3rd: boolean, optional | |
206 | If the input coordinates are 3 dimensional and the 3rd dimension is time. | |
190 | 207 | radians: boolean, optional |
191 | 208 | If True, will expect input data to be in radians and will return radians |
192 | 209 | if the projection is geographic. Default is False (degrees). Ignored for |
216 | 233 | >>> transprojr = Transformer.from_proj('+init=EPSG:4326', {"proj":'geocent', "ellps":'WGS84', "datum":'WGS84'}) |
217 | 234 | >>> for pt in transprojr.itransform([(-2.137, 0.661, -20.531)], radians=True): '{:.3f} {:.3f} {:.3f}'.format(*pt) |
218 | 235 | '-2704214.394 -4254414.478 3894270.731' |
219 | >>> transproj_eq = Transformer.from_proj('+init=EPSG:4326', 4326, skip_equivalent=True) | |
236 | >>> transproj_eq = Transformer.from_proj('+init=EPSG:4326', '+proj=longlat +datum=WGS84 +no_defs +type=crs', skip_equivalent=True) | |
220 | 237 | >>> for pt in transproj_eq.itransform([(-2.137, 0.661)]): '{:.3f} {:.3f}'.format(*pt) |
221 | 238 | '-2.137 0.661' |
222 | 239 | |
229 | 246 | raise ValueError("iterable must contain at least one point") |
230 | 247 | |
231 | 248 | stride = len(fst_pt) |
232 | if stride not in (2, 3): | |
233 | raise ValueError("points can contain up to 3 coordinates") | |
249 | if stride not in (2, 3, 4): | |
250 | raise ValueError("points can contain up to 4 coordinates") | |
251 | ||
252 | if time_3rd and stride != 3: | |
253 | raise ValueError("'time_3rd' is only valid for 3 coordinates.") | |
234 | 254 | |
235 | 255 | # create a coordinate sequence generator etc. x1,y1,z1,x2,y2,z2,.... |
236 | 256 | # chain so the generator returns the first point that was already acquired |
243 | 263 | break |
244 | 264 | |
245 | 265 | self._transformer._transform_sequence( |
246 | stride, buff, switch, radians, errcheck=errcheck | |
266 | stride, | |
267 | buff, | |
268 | switch=switch, | |
269 | time_3rd=time_3rd, | |
270 | radians=radians, | |
271 | errcheck=errcheck, | |
247 | 272 | ) |
248 | 273 | |
249 | 274 | for pt in zip(*([iter(buff)] * stride)): |
251 | 276 | |
252 | 277 | |
253 | 278 | def transform( |
254 | p1, p2, x, y, z=None, radians=False, errcheck=False, skip_equivalent=False | |
279 | p1, p2, x, y, z=None, tt=None, radians=False, errcheck=False, skip_equivalent=False | |
255 | 280 | ): |
256 | 281 | """ |
257 | 282 | x2, y2, z2 = transform(p1, p2, x1, y1, z1) |
333 | 358 | >>> x2, y2 = transform(c1, c2, x1, y1) |
334 | 359 | >>> "%s %s" % (str(x2)[:9],str(y2)[:9]) |
335 | 360 | '1402291.0 5076289.5' |
336 | >>> pj = Proj(init="epsg:4214") | |
337 | >>> pjx, pjy = pj(116.366, 39.867) | |
338 | >>> xr, yr = transform(pj, Proj(4326), pjx, pjy, radians=True, errcheck=True) | |
339 | >>> "%.3f %.3f" % (xr, yr) | |
340 | '0.696 2.031' | |
341 | 361 | >>> xeq, yeq = transform(4326, 4326, 30, 60, skip_equivalent=True) |
342 | 362 | >>> "%.0f %.0f" % (xeq, yeq) |
343 | 363 | '30 60' |
344 | 364 | """ |
345 | 365 | return Transformer.from_proj(p1, p2, skip_equivalent=skip_equivalent).transform( |
346 | x, y, z, radians, errcheck=errcheck | |
366 | xx=x, yy=y, zz=z, tt=tt, radians=radians, errcheck=errcheck | |
347 | 367 | ) |
348 | 368 | |
349 | 369 | |
350 | 370 | def itransform( |
351 | p1, p2, points, switch=False, radians=False, errcheck=False, skip_equivalent=False | |
371 | p1, | |
372 | p2, | |
373 | points, | |
374 | switch=False, | |
375 | time_3rd=False, | |
376 | radians=False, | |
377 | errcheck=False, | |
378 | skip_equivalent=False, | |
352 | 379 | ): |
353 | 380 | """ |
354 | 381 | points2 = itransform(p1, p2, points1) |
395 | 422 | '411050.470 4497928.574' |
396 | 423 | '399060.236 4486978.710' |
397 | 424 | '458553.243 4523045.485' |
398 | >>> pj = Proj(init="epsg:4214") | |
399 | >>> pjx, pjy = pj(116.366, 39.867) | |
400 | >>> for pt in itransform(pj, Proj(4326), [(pjx, pjy)], radians=True, errcheck=True): '{:.3f} {:.3f}'.format(*pt) | |
401 | '0.696 2.031' | |
402 | 425 | >>> for pt in itransform(4326, 4326, [(30, 60)], skip_equivalent=True): '{:.0f} {:.0f}'.format(*pt) |
403 | 426 | '30 60' |
404 | 427 | |
405 | 428 | """ |
406 | 429 | return Transformer.from_proj(p1, p2, skip_equivalent=skip_equivalent).itransform( |
407 | points, switch, radians, errcheck=errcheck | |
430 | points, switch=switch, time_3rd=time_3rd, radians=radians, errcheck=errcheck | |
408 | 431 | ) |
0 | 0 | cython>=0.28.4 |
1 | 1 | mock |
2 | nose2[coverage_plugin]>=0.6.5 | |
3 | pytest | |
4 | cov-core | |
5 | coverage>=4.0 | |
2 | pytest>3.6 | |
3 | pytest-cov | |
6 | 4 | numpy |
13 | 13 | |
14 | 14 | |
15 | 15 | pyproj.datadir.append_data_dir |
16 | --------------------------- | |
16 | ------------------------------ | |
17 | 17 | |
18 | 18 | .. autofunction:: pyproj.datadir.append_data_dir |
0 | 0 | Transformer |
1 | 1 | =========== |
2 | ||
3 | .. warning:: The axis order may be swapped if the source and destination | |
4 | CRS's are defined as having the first coordinate component point in a | |
5 | northerly direction (See PROJ FAQ on | |
6 | `axis order <https://proj4.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent>`_). | |
7 | You can check the axis order with the `pyproj.CRS` class. | |
8 | ||
2 | 9 | |
3 | 10 | pyproj.Transformer |
4 | 11 | ------------------ |
0 | 0 | Change Log |
1 | 1 | ========== |
2 | ||
3 | 2.1.3 | |
4 | ~~~~~ | |
5 | * Added support for time transformations (issue #208) | |
6 | * Fixed projection equivalence testing for transformations (pull #231). | |
7 | * Switch to pytest for testing (pull #230) | |
8 | * Various testing fixes (pull #223, #222, #221, #220) | |
9 | * Convert PROJ error messages from bytes to strings (pull #219) | |
10 | * Fix data dir path separator to be (;) for windows and (:) for linux (pull #234) | |
2 | 11 | |
3 | 12 | 2.1.2 |
4 | 13 | ~~~~~ |
5 | 14 | * Updated to use the CRS definition for Proj instances in transforms (issue #207) |
6 | 15 | * Add option to skip tranformation operation if input and output projections are equivalent |
7 | and always skip if the input and output projections are exact (issue #128) | |
8 | * Update setup.py method for checking PROJ version (pull #211) | |
16 | and always skip if the input and output projections are exact (issue #128) | |
17 | * Update setup.py method for checking PROJ version (pull #211) | |
18 | * Add internal proj error log messages to exceptions (pull #215) | |
9 | 19 | |
10 | 20 | 2.1.1 |
11 | 21 | ~~~~~ |
12 | 22 | * Restore behavior of 1.9.6 when illegal projection transformation requested |
13 | 23 | (return ``inf`` instead of raising an exception, issue #202). kwarg ``errcheck`` |
14 | added to :func:`pyproj.transform` and :func:`pyproj.itransform` | |
24 | added to :func:`~pyproj.transformer.transform` and :func:`~pyproj.transformer.itransform` | |
15 | 25 | (default ``False``). When ``errcheck=True`` an exception is raised. |
16 | 26 | |
17 | 27 | 2.1.0 |
18 | 28 | ~~~~~ |
19 | * Added :class:`pyproj.Transformer` to make repetitive transformations more efficient (issue #187) | |
29 | * Added :class:`~pyproj.transformer.Transformer` to make repetitive transformations more efficient (issue #187) | |
20 | 30 | * Added fix for using local datumgrids with transform (issue #191) |
21 | * Added :class:`pyproj.Transformer.from_pipeline` to support pipeline transformations. | |
31 | * Added :meth:`~pyproj.transformer.Transformer.from_pipeline` to support pipeline transformations. | |
22 | 32 | * Added fix for conversion between radians/degrees for transformations (issues #192 & #195) |
23 | 33 | |
24 | 34 | 2.0.2 |
19 | 19 | y_coords = np.random.randint(200000, 250000) |
20 | 20 | |
21 | 21 | |
22 | Example with `pyproj.transform`: | |
22 | Example with :func:`~pyproj.transformer.transform`: | |
23 | 23 | |
24 | 24 | .. code-block:: python |
25 | 25 | |
27 | 27 | |
28 | 28 | Results: 160 ms ± 3.68 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) |
29 | 29 | |
30 | Example with :class:`pyproj.Transformer`: | |
30 | Example with :class:`~pyproj.transformer.Transformer`: | |
31 | 31 | |
32 | 32 | .. code-block:: python |
33 | 33 |
0 | from pyproj import Proj, transform | |
1 | ||
2 | # illustrates the use of the transform function to | |
3 | # perform coordinate transformations with datum shifts. | |
4 | # | |
5 | # This example is from Roberto Vidmar | |
6 | # | |
7 | # Test point is Trieste, Molo Sartorio | |
8 | # | |
9 | # This data come from the Istituto Geografico Militare (IGM), as well as | |
10 | # the 7 parameters to transform from Gauss-Boaga (our reference frame) | |
11 | # to WGS84 | |
12 | # | |
13 | # WGS84 Lat: 45d38'49.879" (45.647188611) | |
14 | # WGS84 Lon: 13d45'34.397" (13.759554722) | |
15 | # WGS84 z: 52.80; | |
16 | # UTM 33: 403340.97 5055597.17 | |
17 | # GB: 2423346.99 5055619.87 | |
18 | ||
19 | UTM_x = 403340.97 | |
20 | UTM_y = 5055597.17 | |
21 | GB_x = 2423346.99 | |
22 | GB_y = 5055619.87 | |
23 | WGS84_lat = 45.647188611 # Degrees | |
24 | WGS84_lon = 13.759554722 # Degrees | |
25 | UTM_z = WGS84_z = 52.8 # Ellipsoidical height in meters | |
26 | ||
27 | wgs84 = Proj(proj="latlong", datum="WGS84") | |
28 | print("proj4 library version = ", wgs84.proj_version) | |
29 | utm33 = Proj(proj="utm", zone="33") | |
30 | gaussb = Proj( | |
31 | init="epsg:26592", towgs84="-122.74,-34.27,-22.83,-1.884,-3.400,-3.030,-15.62" | |
32 | ) | |
33 | ||
34 | print("\nWGS84-->UTM") | |
35 | print("Trieste, Molo Sartorio WGS84:", WGS84_lon, WGS84_lat, WGS84_z) | |
36 | print("Trieste, Molo Sartorio UTM33 (from IGM):", UTM_x, UTM_y) | |
37 | xutm33, yutm33, zutm33 = transform(wgs84, utm33, WGS84_lon, WGS84_lat, WGS84_z) | |
38 | print("Trieste, Molo Sartorio UTM33 (converted):", xutm33, yutm33, zutm33) | |
39 | print("Difference (meters):", xutm33 - UTM_x, yutm33 - UTM_y, zutm33 - UTM_z) | |
40 | ||
41 | print("\nWGS84-->Gauss-Boaga") | |
42 | print("Trieste, Molo Sartorio Gauss-Boaga (from IGM):", GB_x, GB_y) | |
43 | z = 0 # No ellipsoidical height for Gauss-Boaga | |
44 | xgb, ygb, zgb = transform(wgs84, gaussb, WGS84_lon, WGS84_lat, z) | |
45 | print("Trieste, Molo Sartorio Gauss-Boaga (converted):", xgb, ygb) | |
46 | print("Difference (meters):", xgb - GB_x, ygb - GB_y) | |
47 | ||
48 | print("\nUTM-->WGS84") | |
49 | print("Trieste, Molo Sartorio UTM33 (converted):", xutm33, yutm33, zutm33) | |
50 | back_lon, back_lat, back_z = transform(utm33, wgs84, xutm33, yutm33, zutm33) | |
51 | print("Trieste, Molo Sartorio WGS84 (converted back):", back_lon, back_lat, back_z) | |
52 | print( | |
53 | "Difference (seconds):", | |
54 | (back_lon - WGS84_lon) * 3600, | |
55 | (back_lat - WGS84_lat) * 3600, | |
56 | back_z - WGS84_z, | |
57 | "(m)", | |
58 | ) | |
59 | ||
60 | print("\nGauss-Boaga-->WGS84") | |
61 | print("Trieste, Molo Sartorio Gauss-Boaga (converted):", xgb, ygb, zgb) | |
62 | back_lon, back_lat, back_z = transform(gaussb, wgs84, xgb, ygb, zgb) | |
63 | print("Trieste, Molo Sartorio WGS84 (converted back):", back_lon, back_lat, back_z) | |
64 | print( | |
65 | "Difference (seconds):", | |
66 | (back_lon - WGS84_lon) * 3600, | |
67 | (back_lat - WGS84_lat) * 3600, | |
68 | back_z - WGS84_z, | |
69 | "(m)", | |
70 | ) | |
71 | ||
72 | print("\nUTM (from IGM) --> WGS84") | |
73 | print("Trieste, Molo Sartorio UTM33 (from IGM):", UTM_x, UTM_y) | |
74 | conv_lon, conv_lat, conv_z = transform(utm33, wgs84, UTM_x, UTM_y, UTM_z) | |
75 | print("Trieste, Molo Sartorio WGS84 (converted):", conv_lon, conv_lat, conv_z) | |
76 | print( | |
77 | "Difference (seconds):", | |
78 | (conv_lon - WGS84_lon) * 3600, | |
79 | (conv_lat - WGS84_lat) * 3600, | |
80 | conv_z - WGS84_z, | |
81 | "(m)", | |
82 | ) | |
83 | ||
84 | ||
85 | print("\nGauss-Boaga (from IGM) --> WGS84") | |
86 | print("Trieste, Molo Sartorio Gauss-Boaga (from IGM):", GB_x, GB_y) | |
87 | GB_z = 0 # No ellipsoidical height for Gauss-Boaga | |
88 | conv_lon, conv_lat, conv_z = transform(gaussb, wgs84, GB_x, GB_y, GB_z) | |
89 | print("Trieste, Molo Sartorio WGS84 (converted):", conv_lon, conv_lat) | |
90 | print( | |
91 | "Difference (seconds):", | |
92 | (conv_lon - WGS84_lon) * 3600, | |
93 | (conv_lat - WGS84_lat) * 3600, | |
94 | ) |
0 | from pyproj import Geod | |
1 | import commands, cPickle | |
2 | ||
3 | g = Geod(ellps="clrk66") | |
4 | lat1pt = 42.0 + (15.0 / 60.0) | |
5 | lon1pt = -71.0 - (7.0 / 60.0) | |
6 | lat2pt = 45.0 + (31.0 / 60.0) | |
7 | lon2pt = -123.0 - (41.0 / 60.0) | |
8 | print(lat1pt, lon1pt, lat2pt, lon2pt) | |
9 | print("inverse transform") | |
10 | print("from proj.4 invgeod:") | |
11 | print( | |
12 | commands.getoutput( | |
13 | "echo \"42d15'N 71d07'W 45d31'N 123d41'W\" | geod +ellps=clrk66 -I -f \"%.3f\"" | |
14 | ) | |
15 | ) | |
16 | print("from pyproj.Geod.inv:") | |
17 | az12, az21, dist = g.inv(lon1pt, lat1pt, lon2pt, lat2pt) | |
18 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
19 | print("forward transform") | |
20 | print("from proj.4 geod:") | |
21 | print( | |
22 | commands.getoutput( | |
23 | 'echo "42d15\'N 71d07\'W -66d31\'50.141 4164192.708" | geod +ellps=clrk66 -f "%.3f"' | |
24 | ) | |
25 | ) | |
26 | endlon, endlat, backaz = g.fwd(lon1pt, lat1pt, az12, dist) | |
27 | print("from pyproj.Geod.fwd:") | |
28 | print("%6.3f %6.3f %13.3f" % (endlat, endlon, backaz)) | |
29 | print("intermediate points:") | |
30 | print("from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:") | |
31 | points = "+lon_1=%s +lat_1=%s +lon_2=%s +lat_2=%s" % (lon1pt, lat1pt, lon2pt, lat2pt) | |
32 | print(points) | |
33 | print(commands.getoutput('geod +ellps=clrk66 -f "%.3f" +n_S=5 ' + points)) | |
34 | print("from pyproj.Geod.npts:") | |
35 | npts = 4 | |
36 | lonlats = g.npts(lon1pt, lat1pt, lon2pt, lat2pt, npts) | |
37 | lonprev = lon1pt | |
38 | latprev = lat1pt | |
39 | print(dist / (npts + 1)) | |
40 | print("%6.3f %7.3f" % (lat1pt, lon1pt)) | |
41 | for lon, lat in lonlats: | |
42 | az12, az21, dist = g.inv(lonprev, latprev, lon, lat) | |
43 | print("%6.3f %7.3f %11.3f" % (lat, lon, dist)) | |
44 | latprev = lat | |
45 | lonprev = lon | |
46 | az12, az21, dist = g.inv(lonprev, latprev, lon2pt, lat2pt) | |
47 | print("%6.3f %7.3f %11.3f" % (lat2pt, lon2pt, dist)) | |
48 | ||
49 | # specify the lat/lons of some cities. | |
50 | boston_lat = 42.0 + (15.0 / 60.0) | |
51 | boston_lon = -71.0 - (7.0 / 60.0) | |
52 | portland_lat = 45.0 + (31.0 / 60.0) | |
53 | portland_lon = -123.0 - (41.0 / 60.0) | |
54 | newyork_lat = 40.0 + (47.0 / 60.0) | |
55 | newyork_lon = -73.0 - (58.0 / 60.0) | |
56 | london_lat = 51.0 + (32.0 / 60.0) | |
57 | london_lon = -(5.0 / 60.0) | |
58 | g1 = Geod(ellps="clrk66") | |
59 | cPickle.dump(g1, open("geod1.pickle", "wb"), -1) | |
60 | g2 = Geod(ellps="WGS84") | |
61 | cPickle.dump(g2, open("geod2.pickle", "wb"), -1) | |
62 | az12, az21, dist = g1.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
63 | print("distance between boston and portland, clrk66:") | |
64 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
65 | print("distance between boston and portland, WGS84:") | |
66 | az12, az21, dist = g2.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
67 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
68 | print("testing pickling of Geod instance") | |
69 | g3 = cPickle.load(open("geod1.pickle", "rb")) | |
70 | g4 = cPickle.load(open("geod2.pickle", "rb")) | |
71 | az12, az21, dist = g3.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
72 | print("distance between boston and portland, clrk66 (from pickle):") | |
73 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
74 | az12, az21, dist = g4.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
75 | print("distance between boston and portland, WGS84 (from pickle):") | |
76 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
77 | g3 = Geod("+ellps=clrk66") # proj4 style init string | |
78 | print("inverse transform") | |
79 | print("from proj.4 invgeod:") | |
80 | print( | |
81 | commands.getoutput( | |
82 | "echo \"42d15'N 71d07'W 45d31'N 123d41'W\" | geod +ellps=clrk66 -I -f \"%.3f\"" | |
83 | ) | |
84 | ) | |
85 | print("from pyproj.Geod.inv:") | |
86 | az12, az21, dist = g3.inv(lon1pt, lat1pt, lon2pt, lat2pt) | |
87 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) |
0 | from pyproj import Proj | |
1 | import time, cPickle, array | |
2 | import numpy as N | |
3 | ||
4 | params = {} | |
5 | params["proj"] = "lcc" | |
6 | params["R"] = 6371200 | |
7 | params["lat_1"] = 50 | |
8 | params["lat_2"] = 50 | |
9 | params["lon_0"] = -107 | |
10 | nx = 349 | |
11 | ny = 277 | |
12 | dx = 32463.41 | |
13 | dy = dx | |
14 | # can either use a dict | |
15 | # awips221 = Proj(params) | |
16 | # or keyword args | |
17 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
18 | print("proj4 library version = ", awips221.proj_version) | |
19 | # AWIPS grid 221 parameters | |
20 | # (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html) | |
21 | llcrnrx, llcrnry = awips221(-145.5, 1.0) | |
22 | params["x_0"] = -llcrnrx | |
23 | params["y_0"] = -llcrnry | |
24 | awips221 = Proj(params) | |
25 | llcrnrx, llcrnry = awips221(-145.5, 1.0) | |
26 | # find 4 lon/lat crnrs of AWIPS grid 221. | |
27 | llcrnrx = 0.0 | |
28 | llcrnry = 0.0 | |
29 | lrcrnrx = dx * (nx - 1) | |
30 | lrcrnry = 0.0 | |
31 | ulcrnrx = 0.0 | |
32 | ulcrnry = dy * (ny - 1) | |
33 | urcrnrx = dx * (nx - 1) | |
34 | urcrnry = dy * (ny - 1) | |
35 | llcrnrlon, llcrnrlat = awips221(llcrnrx, llcrnry, inverse=True) | |
36 | lrcrnrlon, lrcrnrlat = awips221(lrcrnrx, lrcrnry, inverse=True) | |
37 | urcrnrlon, urcrnrlat = awips221(urcrnrx, urcrnry, inverse=True) | |
38 | ulcrnrlon, ulcrnrlat = awips221(ulcrnrx, ulcrnry, inverse=True) | |
39 | print("4 crnrs of AWIPS grid 221:") | |
40 | print(llcrnrlon, llcrnrlat) | |
41 | print(lrcrnrlon, lrcrnrlat) | |
42 | print(urcrnrlon, urcrnrlat) | |
43 | print(ulcrnrlon, ulcrnrlat) | |
44 | print("from GRIB docs") | |
45 | print("(see http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)") | |
46 | print(" -145.5 1.0") | |
47 | print(" -68.318 0.897") | |
48 | print(" -2.566 46.352") | |
49 | print(" 148.639 46.635") | |
50 | # compute lons and lats for the whole AWIPS grid 221 (377x249). | |
51 | dx = (urcrnrx - llcrnrx) / (nx - 1) | |
52 | dy = (urcrnry - llcrnry) / (ny - 1) | |
53 | x = llcrnrx + dx * N.indices((ny, nx), "f")[1, :, :] | |
54 | y = llcrnry + dy * N.indices((ny, nx), "f")[0, :, :] | |
55 | t1 = time.clock() | |
56 | lons, lats = awips221(N.ravel(x).tolist(), N.ravel(y).tolist(), inverse=True) | |
57 | t2 = time.clock() | |
58 | print("data in lists:") | |
59 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
60 | print("max/min lons") | |
61 | print(min(lons), max(lons)) | |
62 | print("max/min lats") | |
63 | print(min(lats), max(lats)) | |
64 | print("took", t2 - t1, "secs") | |
65 | xa = array.array("f", N.ravel(x).tolist()) | |
66 | ya = array.array("f", N.ravel(y).tolist()) | |
67 | t1 = time.clock() | |
68 | lons, lats = awips221(xa, ya, inverse=True) | |
69 | t2 = time.clock() | |
70 | print("data in python arrays:") | |
71 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
72 | print("max/min lons") | |
73 | print(min(lons), max(lons)) | |
74 | print("max/min lats") | |
75 | print(min(lats), max(lats)) | |
76 | print("took", t2 - t1, "secs") | |
77 | t1 = time.clock() | |
78 | lons, lats = awips221(x, y, inverse=True) | |
79 | t2 = time.clock() | |
80 | print("data in a numpy array:") | |
81 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
82 | print("max/min lons") | |
83 | print(N.minimum.reduce(N.ravel(lons)), N.maximum.reduce(N.ravel(lons))) | |
84 | print("max/min lats") | |
85 | print(N.minimum.reduce(N.ravel(lats)), N.maximum.reduce(N.ravel(lats))) | |
86 | print("took", t2 - t1, "secs") | |
87 | # reverse transformation. | |
88 | t1 = time.clock() | |
89 | xx, yy = awips221(lons, lats) | |
90 | t2 = time.clock() | |
91 | print("max abs error for x") | |
92 | print(N.maximum.reduce(N.fabs(N.ravel(x - xx)))) | |
93 | print("max abs error for y") | |
94 | print(N.maximum.reduce(N.fabs(N.ravel(x - xx)))) | |
95 | print("took", t2 - t1, "secs") | |
96 | cPickle.dump(awips221, open("test.pickle", "wb"), -1) | |
97 | print("compare output with sample.out") | |
98 | print("now run test2.py to test pickling") | |
0 | # -*- coding: utf-8 -*- | |
1 | """Rewrite part of test.py in pyproj in the form of unittests.""" | |
2 | ||
3 | import math | |
4 | import sys | |
5 | import unittest | |
6 | ||
7 | from pyproj import Geod, Proj, pj_ellps, pj_list, transform | |
8 | from pyproj.crs import CRSError | |
9 | ||
10 | ||
11 | class BasicTest(unittest.TestCase): | |
12 | def testProj4Version(self): | |
13 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
14 | # self.assertEqual(awips221.proj_version, 4.9) | |
15 | ||
16 | def testInitWithBackupString4(self): | |
17 | # this fails unless backup of to_string(4) is used | |
18 | pj = Proj( | |
19 | "+proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@null +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs" | |
20 | ) | |
21 | assert pj.crs.is_valid | |
22 | ||
23 | def testProjAwips221(self): | |
24 | # AWIPS is Advanced Weather Interactive Processing System | |
25 | params = {"proj": "lcc", "R": 6371200, "lat_1": 50, "lat_2": 50, "lon_0": -107} | |
26 | nx = 349 | |
27 | ny = 277 | |
28 | awips221 = Proj( | |
29 | proj=params["proj"], | |
30 | R=params["R"], | |
31 | lat_1=params["lat_1"], | |
32 | lat_2=params["lat_2"], | |
33 | lon_0=params["lon_0"], | |
34 | preserve_units=False, | |
35 | ) | |
36 | awips221_from_dict = Proj(params, preserve_units=False) | |
37 | ||
38 | items = sorted([val for val in awips221.crs.srs.split() if val]) | |
39 | items_dict = sorted([val for val in awips221_from_dict.crs.srs.split() if val]) | |
40 | self.assertEqual(items, items_dict) | |
41 | ||
42 | expected = sorted( | |
43 | [ | |
44 | "+proj=lcc", | |
45 | "+R=6371200", | |
46 | "+lat_1=50", | |
47 | "+lat_2=50", | |
48 | "+lon_0=-107", | |
49 | "+type=crs", | |
50 | ] | |
51 | ) | |
52 | self.assertEqual(items, expected) | |
53 | ||
54 | point = awips221(-145.5, 1.0) | |
55 | x, y = -5632642.22547495, 1636571.4883145525 | |
56 | self.assertAlmostEqual(point[0], x) | |
57 | self.assertAlmostEqual(point[1], y) | |
58 | ||
59 | pairs = [ | |
60 | [(-45, 45), (4351601.20766915, 7606948.029327129)], | |
61 | [(45, 45), (5285389.07739382, 14223336.17467613)], | |
62 | [(45, -45), (20394982.466924712, 21736546.456803113)], | |
63 | [(-45, -45), (16791730.756976362, -3794425.4816524936)], | |
64 | ] | |
65 | for point_geog, expected in pairs: | |
66 | point = awips221(*point_geog) | |
67 | self.assertAlmostEqual(point[0], expected[0]) | |
68 | self.assertAlmostEqual(point[1], expected[1]) | |
69 | point_geog2 = awips221(*point, inverse=True) | |
70 | self.assertAlmostEqual(point_geog[0], point_geog2[0]) | |
71 | self.assertAlmostEqual(point_geog[1], point_geog2[1]) | |
72 | ||
73 | def test_from_dict_with_bool(self): | |
74 | # issue #183 | |
75 | p_d = {'proj': 'omerc', | |
76 | 'lat_2': 80.27942, | |
77 | 'lat_0': 62.87671, | |
78 | 'lat_1': 42.751232, | |
79 | 'ellps': 'WGS84', | |
80 | 'no_rot': True, | |
81 | 'lon_1': 33.793186, | |
82 | 'lon_2': -18.374414} | |
83 | p=Proj(p_d) | |
84 | self.assertTrue('+no_rot' in p.srs.split()) | |
85 | p_d = {'proj': 'omerc', | |
86 | 'lat_2': 80.27942, | |
87 | 'lat_0': 62.87671, | |
88 | 'lat_1': 42.751232, | |
89 | 'ellps': 'WGS84', | |
90 | 'no_rot': False, | |
91 | 'lon_1': 33.793186, | |
92 | 'lon_2': -18.374414} | |
93 | p=Proj(p_d) | |
94 | self.assertFalse('+no_rot' in p.srs.split()) | |
95 | ||
96 | ||
97 | ||
98 | class InverseHammerTest(unittest.TestCase): | |
99 | # This is a unit test of the inverse of the hammer projection, which | |
100 | # was added in the 4.9.3 version of PROJ (then PROJ.4). | |
101 | @classmethod | |
102 | def setUpClass(self): | |
103 | self.p = Proj(proj="hammer") # hammer proj | |
104 | self.x, self.y = self.p(-30, 40) | |
105 | ||
106 | def test_forward(self): | |
107 | self.assertAlmostEqual(self.x, -2711575.083, places=3) | |
108 | self.assertAlmostEqual(self.y, 4395506.619, places=3) | |
109 | ||
110 | def test_inverse(self): | |
111 | lon, lat = self.p(self.x, self.y, inverse=True) | |
112 | self.assertAlmostEqual(lon, -30.0, places=3) | |
113 | self.assertAlmostEqual(lat, 40.0, places=3) | |
114 | ||
115 | ||
116 | class TypeError_Transform_Issue8_Test(unittest.TestCase): | |
117 | # Test for "Segmentation fault on pyproj.transform #8" | |
118 | # https://github.com/jswhit/pyproj/issues/8 | |
119 | ||
120 | def setUp(self): | |
121 | self.p = Proj(init="epsg:4269") | |
122 | ||
123 | def test_tranform_none_1st_parmeter(self): | |
124 | # test should raise Type error if projections are not of Proj classes | |
125 | # version 1.9.4 produced AttributeError, now should raise TypeError | |
126 | with self.assertRaises(CRSError): | |
127 | transform(None, self.p, -74, 39) | |
128 | ||
129 | def test_tranform_none_2nd_parmeter(self): | |
130 | # test should raise Type error if projections are not of Proj classes | |
131 | # version 1.9.4 has a Segmentation Fault, now should raise TypeError | |
132 | with self.assertRaises(CRSError): | |
133 | transform(self.p, None, -74, 39) | |
134 | ||
135 | ||
136 | class Geod_NoDefs_Issue22_Test(unittest.TestCase): | |
137 | # Test for Issue #22, Geod with "+no_defs" in initstring | |
138 | # Before PR #23 merged 2015-10-07, having +no_defs in the initstring would result in a ValueError | |
139 | def test_geod_nodefs(self): | |
140 | Geod("+a=6378137 +b=6378137 +no_defs") | |
141 | ||
142 | ||
143 | class ProjLatLongTypeErrorTest(unittest.TestCase): | |
144 | # .latlong() using in transform raised a TypeError in release 1.9.5.1 | |
145 | # reported in issue #53, resolved in #73. | |
146 | def test_latlong_typeerror(self): | |
147 | p = Proj("+proj=stere +lon_0=-39 +lat_0=90 +lat_ts=71.0 +ellps=WGS84") | |
148 | self.assertTrue(isinstance(p, Proj)) | |
149 | # if not patched this line raises a "TypeError: p2 must be a Proj class" | |
150 | lon, lat = transform(p, p.to_latlong(), 200000, 400000) | |
151 | ||
152 | ||
153 | @unittest.skipIf(sys.version_info < (3, 4), "Python 3.4 or newer required for subTest()") | |
154 | class ForwardInverseTest(unittest.TestCase): | |
155 | def test_fwd_inv(self): | |
156 | for pj in pj_list.keys(): | |
157 | with self.subTest(pj=pj): | |
158 | try: | |
159 | p = Proj(proj=pj) | |
160 | x, y = p(-30, 40) | |
161 | # note, for proj 4.9.2 or before the inverse projection may be missing | |
162 | # and pyproj 1.9.5.1 or before does not test for this and will | |
163 | # give a segmentation fault at this point: | |
164 | lon, lat = p(x, y, inverse=True) | |
165 | except RuntimeError: | |
166 | pass | |
167 | ||
168 | ||
169 | # Tests for shared memory between Geod objects | |
170 | class GeodSharedMemoryBugTestIssue64(unittest.TestCase): | |
171 | def setUp(self): | |
172 | self.g = Geod(ellps="clrk66") | |
173 | self.ga = self.g.a | |
174 | self.mercury = Geod(a=2439700) # Mercury 2000 ellipsoid | |
175 | # Mercury is much smaller than earth. | |
176 | ||
177 | def test_not_shared_memory(self): | |
178 | self.assertEqual(self.ga, self.g.a) | |
179 | # mecury must have a different major axis from earth | |
180 | self.assertNotEqual(self.g.a, self.mercury.a) | |
181 | self.assertNotEqual(self.g.b, self.mercury.b) | |
182 | self.assertNotEqual(self.g.sphere, self.mercury.sphere) | |
183 | self.assertNotEqual(self.g.f, self.mercury.f) | |
184 | self.assertNotEqual(self.g.es, self.mercury.es) | |
185 | ||
186 | # initstrings were not shared in issue #64 | |
187 | self.assertNotEqual(self.g.initstring, self.mercury.initstring) | |
188 | ||
189 | def test_distances(self): | |
190 | # note calculated distance was not an issue with #64, but it still a shared memory test | |
191 | boston_lat = 42.0 + (15.0 / 60.0) | |
192 | boston_lon = -71.0 - (7.0 / 60.0) | |
193 | portland_lat = 45.0 + (31.0 / 60.0) | |
194 | portland_lon = -123.0 - (41.0 / 60.0) | |
195 | ||
196 | az12, az21, dist_g = self.g.inv( | |
197 | boston_lon, boston_lat, portland_lon, portland_lat | |
198 | ) | |
199 | ||
200 | az12, az21, dist_mercury = self.mercury.inv( | |
201 | boston_lon, boston_lat, portland_lon, portland_lat | |
202 | ) | |
203 | self.assertLess(dist_mercury, dist_g) | |
204 | ||
205 | ||
206 | class ReprTests(unittest.TestCase): | |
207 | # test __repr__ for Proj object | |
208 | def test_repr(self): | |
209 | p = Proj(proj="latlong", preserve_units=True) | |
210 | expected = "Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True)" | |
211 | self.assertEqual(repr(p), expected) | |
212 | ||
213 | # test __repr__ for Geod object | |
214 | def test_sphere(self): | |
215 | # ellipse is Venus 2000 (IAU2000:29900), which is a sphere | |
216 | g = Geod("+a=6051800 +b=6051800") | |
217 | self.assertEqual(repr(g), "Geod('+a=6051800 +f=0')") | |
218 | ||
219 | # test __repr__ for Geod object | |
220 | def test_ellps_name_round_trip(self): | |
221 | # this could be done in a parameter fashion | |
222 | for ellps_name in pj_ellps: | |
223 | # skip tests, these ellipses NWL9D and WGS66 are the same | |
224 | if ellps_name in ("NWL9D", "WGS66"): | |
225 | continue | |
226 | p = Geod(ellps=ellps_name) | |
227 | expected = "Geod(ellps='{0}')".format(ellps_name) | |
228 | self.assertEqual(repr(p), expected) | |
229 | ||
230 | ||
231 | class TestRadians(unittest.TestCase): | |
232 | """Tests issue #84""" | |
233 | ||
234 | def setUp(self): | |
235 | self.g = Geod(ellps="clrk66") | |
236 | self.boston_d = (-71.0 - (7.0 / 60.0), 42.0 + (15.0 / 60.0)) | |
237 | self.boston_r = (math.radians(self.boston_d[0]), math.radians(self.boston_d[1])) | |
238 | self.portland_d = (-123.0 - (41.0 / 60.0), 45.0 + (31.0 / 60.0)) | |
239 | self.portland_r = ( | |
240 | math.radians(self.portland_d[0]), | |
241 | math.radians(self.portland_d[1]), | |
242 | ) | |
243 | ||
244 | def test_inv_radians(self): | |
245 | ||
246 | # Get bearings and distance from Boston to Portland in degrees | |
247 | az12_d, az21_d, dist_d = self.g.inv( | |
248 | self.boston_d[0], | |
249 | self.boston_d[1], | |
250 | self.portland_d[0], | |
251 | self.portland_d[1], | |
252 | radians=False, | |
253 | ) | |
254 | ||
255 | # Get bearings and distance from Boston to Portland in radians | |
256 | az12_r, az21_r, dist_r = self.g.inv( | |
257 | self.boston_r[0], | |
258 | self.boston_r[1], | |
259 | self.portland_r[0], | |
260 | self.portland_r[1], | |
261 | radians=True, | |
262 | ) | |
263 | ||
264 | # Check they are equal | |
265 | self.assertAlmostEqual(az12_d, math.degrees(az12_r)) | |
266 | self.assertAlmostEqual(az21_d, math.degrees(az21_r)) | |
267 | self.assertAlmostEqual(dist_d, dist_r) | |
268 | ||
269 | def test_fwd_radians(self): | |
270 | # Get bearing and distance to Portland | |
271 | az12_d, az21_d, dist = self.g.inv( | |
272 | self.boston_d[0], | |
273 | self.boston_d[1], | |
274 | self.portland_d[0], | |
275 | self.portland_d[1], | |
276 | radians=False, | |
277 | ) | |
278 | ||
279 | # Calculate Portland's lon/lat from bearing and distance in degrees | |
280 | endlon_d, endlat_d, backaz_d = self.g.fwd( | |
281 | self.boston_d[0], self.boston_d[1], az12_d, dist, radians=False | |
282 | ) | |
283 | ||
284 | # Calculate Portland's lon/lat from bearing and distance in radians | |
285 | endlon_r, endlat_r, backaz_r = self.g.fwd( | |
286 | self.boston_r[0], self.boston_r[1], math.radians(az12_d), dist, radians=True | |
287 | ) | |
288 | ||
289 | # Check they are equal | |
290 | self.assertAlmostEqual(endlon_d, math.degrees(endlon_r)) | |
291 | self.assertAlmostEqual(endlat_d, math.degrees(endlat_r)) | |
292 | self.assertAlmostEqual(backaz_d, math.degrees(backaz_r)) | |
293 | ||
294 | # Check to make sure we're back in Portland | |
295 | self.assertAlmostEqual(endlon_d, self.portland_d[0]) | |
296 | self.assertAlmostEqual(endlat_d, self.portland_d[1]) | |
297 | ||
298 | def test_npts_radians(self): | |
299 | # Calculate 10 points between Boston and Portland in degrees | |
300 | points_d = self.g.npts( | |
301 | lon1=self.boston_d[0], | |
302 | lat1=self.boston_d[1], | |
303 | lon2=self.portland_d[0], | |
304 | lat2=self.portland_d[1], | |
305 | npts=10, | |
306 | radians=False, | |
307 | ) | |
308 | ||
309 | # Calculate 10 points between Boston and Portland in radians | |
310 | points_r = self.g.npts( | |
311 | lon1=self.boston_r[0], | |
312 | lat1=self.boston_r[1], | |
313 | lon2=self.portland_r[0], | |
314 | lat2=self.portland_r[1], | |
315 | npts=10, | |
316 | radians=True, | |
317 | ) | |
318 | ||
319 | # Check they are equal | |
320 | for index, dpoint in enumerate(points_d): | |
321 | self.assertAlmostEqual(dpoint[0], math.degrees(points_r[index][0])) | |
322 | self.assertAlmostEqual(dpoint[1], math.degrees(points_r[index][1])) | |
323 | ||
324 | ||
325 | class Geod_NaN_Issue112_Test(unittest.TestCase): | |
326 | # Test for Issue #112; Geod should silently propagate NaNs in input | |
327 | # to the output. | |
328 | def test_geod_nans(self): | |
329 | g = Geod(ellps="clrk66") | |
330 | (azi1, azi2, s12) = g.inv(43, 10, float("nan"), 20) | |
331 | self.assertTrue(azi1 != azi1) | |
332 | self.assertTrue(azi2 != azi2) | |
333 | self.assertTrue(s12 != s12) | |
334 | (azi1, azi2, s12) = g.inv(43, 10, 53, float("nan")) | |
335 | self.assertTrue(azi1 != azi1) | |
336 | self.assertTrue(azi2 != azi2) | |
337 | self.assertTrue(s12 != s12) | |
338 | # Illegal latitude is treated as NaN | |
339 | (azi1, azi2, s12) = g.inv(43, 10, 53, 91) | |
340 | self.assertTrue(azi1 != azi1) | |
341 | self.assertTrue(azi2 != azi2) | |
342 | self.assertTrue(s12 != s12) | |
343 | (lon2, lat2, azi2) = g.fwd(43, 10, float("nan"), 1e6) | |
344 | self.assertTrue(lon2 != lon2) | |
345 | self.assertTrue(lat2 != lat2) | |
346 | self.assertTrue(azi2 != azi2) | |
347 | (lon2, lat2, azi2) = g.fwd(43, 10, 20, float("nan")) | |
348 | self.assertTrue(lon2 != lon2) | |
349 | self.assertTrue(lat2 != lat2) | |
350 | self.assertTrue(azi2 != azi2) | |
351 | (lon2, lat2, azi2) = g.fwd(43, float("nan"), 20, 1e6) | |
352 | self.assertTrue(lon2 != lon2) | |
353 | self.assertTrue(lat2 != lat2) | |
354 | self.assertTrue(azi2 != azi2) | |
355 | # Illegal latitude is treated as NaN | |
356 | (lon2, lat2, azi2) = g.fwd(43, 91, 20, 1e6) | |
357 | self.assertTrue(lon2 != lon2) | |
358 | self.assertTrue(lat2 != lat2) | |
359 | self.assertTrue(azi2 != azi2) | |
360 | # Only lon2 is NaN | |
361 | (lon2, lat2, azi2) = g.fwd(float("nan"), 10, 20, 1e6) | |
362 | self.assertTrue(lon2 != lon2) | |
363 | self.assertTrue(lat2 == lat2) | |
364 | self.assertTrue(azi2 == azi2) | |
365 | ||
366 | ||
367 | def test_proj_equals(): | |
368 | assert Proj(4326) == Proj("epsg:4326") | |
369 | assert Proj(4326) != Proj("epsg:3857") | |
370 | assert Proj(4326) == Proj(Proj("epsg:4326").crs.to_proj4()) | |
371 | ||
372 | ||
373 | if __name__ == "__main__": | |
374 | unittest.main() |
0 | """run test.py first!""" | |
1 | from pyproj import Proj | |
2 | import time, cPickle | |
3 | import numpy as N | |
4 | ||
5 | nx = 349 | |
6 | ny = 277 | |
7 | dx = 32463.41 | |
8 | dy = dx | |
9 | print("do it again, from pickled instance ...") | |
10 | # find 4 lon/lat crnrs of AWIPS grid 221. | |
11 | llcrnrx = 0.0 | |
12 | llcrnry = 0.0 | |
13 | lrcrnrx = dx * (nx - 1) | |
14 | lrcrnry = 0.0 | |
15 | ulcrnrx = 0.0 | |
16 | ulcrnry = dy * (ny - 1) | |
17 | urcrnrx = dx * (nx - 1) | |
18 | urcrnry = dy * (ny - 1) | |
19 | dx = (urcrnrx - llcrnrx) / (nx - 1) | |
20 | dy = (urcrnry - llcrnry) / (ny - 1) | |
21 | x = llcrnrx + dx * N.indices((ny, nx), "f")[1, :, :] | |
22 | y = llcrnry + dy * N.indices((ny, nx), "f")[0, :, :] | |
23 | awips221 = cPickle.load(open("test.pickle", "rb")) | |
24 | t1 = time.clock() | |
25 | lons, lats = awips221(x, y, inverse=True, radians=True) | |
26 | t2 = time.clock() | |
27 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
28 | print("max/min lons in radians") | |
29 | print(N.minimum.reduce(N.ravel(lons)), N.maximum.reduce(N.ravel(lons))) | |
30 | print("max/min lats in radians") | |
31 | print(N.minimum.reduce(N.ravel(lats)), N.maximum.reduce(N.ravel(lats))) | |
32 | print("took", t2 - t1, "secs") | |
33 | # reverse transformation. | |
34 | t1 = time.clock() | |
35 | xx, yy = awips221(lons, lats, radians=True) | |
36 | t2 = time.clock() | |
37 | print("max abs error for x") | |
38 | print(N.maximum.reduce(N.fabs(N.ravel(x - xx)))) | |
39 | print("max abs error for y") | |
40 | print(N.maximum.reduce(N.fabs(N.ravel(y - yy)))) | |
41 | print("took", t2 - t1, "secs") |
0 | import array | |
1 | import time | |
2 | ||
3 | import numpy | |
4 | from numpy.testing import assert_allclose | |
5 | ||
6 | from pyproj import Proj | |
7 | ||
8 | ||
9 | def test_awips221(): | |
10 | params = {} | |
11 | params["proj"] = "lcc" | |
12 | params["R"] = 6371200 | |
13 | params["lat_1"] = 50 | |
14 | params["lat_2"] = 50 | |
15 | params["lon_0"] = -107 | |
16 | nx = 349 | |
17 | ny = 277 | |
18 | dx = 32463.41 | |
19 | dy = dx | |
20 | # can either use a dict | |
21 | # awips221 = Proj(params) | |
22 | # or keyword args | |
23 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
24 | print("proj4 library version = ", awips221.proj_version) | |
25 | # AWIPS grid 221 parameters | |
26 | # (from http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html) | |
27 | llcrnrx, llcrnry = awips221(-145.5, 1.0) | |
28 | params["x_0"] = -llcrnrx | |
29 | params["y_0"] = -llcrnry | |
30 | awips221 = Proj(params) | |
31 | llcrnrx, llcrnry = awips221(-145.5, 1.0) | |
32 | # find 4 lon/lat crnrs of AWIPS grid 221. | |
33 | llcrnrx = 0.0 | |
34 | llcrnry = 0.0 | |
35 | lrcrnrx = dx * (nx - 1) | |
36 | lrcrnry = 0.0 | |
37 | ulcrnrx = 0.0 | |
38 | ulcrnry = dy * (ny - 1) | |
39 | urcrnrx = dx * (nx - 1) | |
40 | urcrnry = dy * (ny - 1) | |
41 | llcrnrlon, llcrnrlat = awips221(llcrnrx, llcrnry, inverse=True) | |
42 | lrcrnrlon, lrcrnrlat = awips221(lrcrnrx, lrcrnry, inverse=True) | |
43 | urcrnrlon, urcrnrlat = awips221(urcrnrx, urcrnry, inverse=True) | |
44 | ulcrnrlon, ulcrnrlat = awips221(ulcrnrx, ulcrnry, inverse=True) | |
45 | print("4 crnrs of AWIPS grid 221:") | |
46 | print(llcrnrlon, llcrnrlat) | |
47 | print(lrcrnrlon, lrcrnrlat) | |
48 | print(urcrnrlon, urcrnrlat) | |
49 | print(ulcrnrlon, ulcrnrlat) | |
50 | print("from GRIB docs") | |
51 | print("(see http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html)") | |
52 | print(" -145.5 1.0") | |
53 | print(" -68.318 0.897") | |
54 | print(" -2.566 46.352") | |
55 | print(" 148.639 46.635") | |
56 | # compute lons and lats for the whole AWIPS grid 221 (377x249). | |
57 | dx = (urcrnrx - llcrnrx) / (nx - 1) | |
58 | dy = (urcrnry - llcrnry) / (ny - 1) | |
59 | x = llcrnrx + dx * numpy.indices((ny, nx), "f")[1, :, :] | |
60 | y = llcrnry + dy * numpy.indices((ny, nx), "f")[0, :, :] | |
61 | t1 = time.clock() | |
62 | lons, lats = awips221( | |
63 | numpy.ravel(x).tolist(), numpy.ravel(y).tolist(), inverse=True | |
64 | ) | |
65 | t2 = time.clock() | |
66 | print("data in lists:") | |
67 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
68 | print("max/min lons") | |
69 | print(min(lons), max(lons)) | |
70 | print("max/min lats") | |
71 | print(min(lats), max(lats)) | |
72 | print("took", t2 - t1, "secs") | |
73 | xa = array.array("f", numpy.ravel(x).tolist()) | |
74 | ya = array.array("f", numpy.ravel(y).tolist()) | |
75 | t1 = time.clock() | |
76 | lons, lats = awips221(xa, ya, inverse=True) | |
77 | t2 = time.clock() | |
78 | print("data in python arrays:") | |
79 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
80 | print("max/min lons") | |
81 | print(min(lons), max(lons)) | |
82 | print("max/min lats") | |
83 | print(min(lats), max(lats)) | |
84 | print("took", t2 - t1, "secs") | |
85 | t1 = time.clock() | |
86 | lons, lats = awips221(x, y, inverse=True) | |
87 | t2 = time.clock() | |
88 | print("data in a numpy array:") | |
89 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
90 | print("max/min lons") | |
91 | print( | |
92 | numpy.minimum.reduce(numpy.ravel(lons)), numpy.maximum.reduce(numpy.ravel(lons)) | |
93 | ) | |
94 | print("max/min lats") | |
95 | print( | |
96 | numpy.minimum.reduce(numpy.ravel(lats)), numpy.maximum.reduce(numpy.ravel(lats)) | |
97 | ) | |
98 | print("took", t2 - t1, "secs") | |
99 | # reverse transformation. | |
100 | t1 = time.clock() | |
101 | xx, yy = awips221(lons, lats) | |
102 | t2 = time.clock() | |
103 | print("max abs error for x") | |
104 | max_abs_err_x = numpy.maximum.reduce(numpy.fabs(numpy.ravel(x - xx))) | |
105 | print(max_abs_err_x) | |
106 | assert_allclose(max_abs_err_x, 0, atol=1e-4) | |
107 | print("max abs error for y") | |
108 | max_abs_err_y = numpy.maximum.reduce(numpy.fabs(numpy.ravel(y - yy))) | |
109 | print(max_abs_err_y) | |
110 | assert_allclose(max_abs_err_y, 0, atol=1e-4) | |
111 | print("took", t2 - t1, "secs") | |
112 | print("compare output with sample.out") |
0 | import pytest | |
1 | ||
2 | from pyproj import CRS | |
3 | from pyproj.exceptions import CRSError | |
4 | ||
5 | ||
6 | def test_from_proj4_json(): | |
7 | json_str = '{"proj": "longlat", "ellps": "WGS84", "datum": "WGS84"}' | |
8 | proj = CRS.from_string(json_str) | |
9 | assert proj.to_proj4(4) == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
10 | assert proj.to_proj4(5) == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
11 | # Test with invalid JSON code | |
12 | with pytest.raises(CRSError): | |
13 | assert CRS.from_string("{foo: bar}") | |
14 | ||
15 | ||
16 | def test_from_epsg(): | |
17 | proj = CRS.from_epsg(4326) | |
18 | assert proj.to_epsg() == 4326 | |
19 | ||
20 | # Test with invalid EPSG code | |
21 | with pytest.raises(CRSError): | |
22 | assert CRS.from_epsg(0) | |
23 | ||
24 | ||
25 | def test_from_epsg_string(): | |
26 | proj = CRS.from_string("epsg:4326") | |
27 | assert proj.to_epsg() == 4326 | |
28 | ||
29 | # Test with invalid EPSG code | |
30 | with pytest.raises(ValueError): | |
31 | assert CRS.from_string("epsg:xyz") | |
32 | ||
33 | ||
34 | def test_from_string(): | |
35 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
36 | assert wgs84_crs.to_proj4() == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
37 | # Make sure this doesn't get handled using the from_epsg() even though 'epsg' is in the string | |
38 | epsg_init_crs = CRS.from_string("+init=epsg:26911 +units=m +no_defs=True") | |
39 | assert ( | |
40 | epsg_init_crs.to_proj4() | |
41 | == "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +type=crs" | |
42 | ) | |
43 | ||
44 | ||
45 | def test_bare_parameters(): | |
46 | """ Make sure that bare parameters (e.g., no_defs) are handled properly, | |
47 | even if they come in with key=True. This covers interaction with pyproj, | |
48 | which makes presents bare parameters as key=<bool>.""" | |
49 | ||
50 | # Example produced by pyproj | |
51 | proj = CRS.from_string( | |
52 | "+proj=lcc +lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
53 | ) | |
54 | assert "+no_defs" in proj.to_proj4(4) | |
55 | ||
56 | # TODO: THIS DOES NOT WORK | |
57 | proj = CRS.from_string( | |
58 | "+lon_0=-95 +ellps=GRS80 +proj=lcc +y_0=0 +no_defs=False +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
59 | ) | |
60 | # assert "+no_defs" not in proj.to_proj4(4) | |
61 | ||
62 | ||
63 | def test_is_geographic(): | |
64 | assert CRS({"init": "EPSG:4326"}).is_geographic is True | |
65 | assert CRS({"init": "EPSG:3857"}).is_geographic is False | |
66 | ||
67 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
68 | assert wgs84_crs.is_geographic is True | |
69 | ||
70 | nad27_crs = CRS.from_string("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs") | |
71 | assert nad27_crs.is_geographic is True | |
72 | ||
73 | lcc_crs = CRS.from_string( | |
74 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
75 | ) | |
76 | assert lcc_crs.is_geographic is False | |
77 | ||
78 | ||
79 | def test_is_projected(): | |
80 | assert CRS({"init": "EPSG:3857"}).is_projected is True | |
81 | ||
82 | lcc_crs = CRS.from_string( | |
83 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
84 | ) | |
85 | assert CRS.from_user_input(lcc_crs).is_projected is True | |
86 | ||
87 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
88 | assert CRS.from_user_input(wgs84_crs).is_projected is False | |
89 | ||
90 | ||
91 | def test_is_same_crs(): | |
92 | crs1 = CRS({"init": "EPSG:4326"}) | |
93 | crs2 = CRS({"init": "EPSG:3857"}) | |
94 | ||
95 | assert crs1 == crs1 | |
96 | assert crs1 != crs2 | |
97 | ||
98 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84") | |
99 | assert crs1 == wgs84_crs | |
100 | ||
101 | # Make sure that same projection with different parameter are not equal | |
102 | lcc_crs1 = CRS.from_string( | |
103 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
104 | ) | |
105 | lcc_crs2 = CRS.from_string( | |
106 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=45 +lat_0=0" | |
107 | ) | |
108 | assert lcc_crs1 != lcc_crs2 | |
109 | ||
110 | ||
111 | def test_to_proj4(): | |
112 | assert ( | |
113 | CRS({"init": "EPSG:4326"}).to_proj4(4) | |
114 | == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
115 | ) | |
116 | ||
117 | ||
118 | def test_is_valid_false(): | |
119 | with pytest.raises(CRSError): | |
120 | CRS(init="EPSG:432600").is_valid | |
121 | ||
122 | ||
123 | def test_is_valid(): | |
124 | assert CRS(init="EPSG:4326").is_valid | |
125 | ||
126 | ||
127 | def test_empty_json(): | |
128 | with pytest.raises(CRSError): | |
129 | CRS.from_string("{}") | |
130 | with pytest.raises(CRSError): | |
131 | CRS.from_string("[]") | |
132 | with pytest.raises(CRSError): | |
133 | CRS.from_string("") | |
134 | ||
135 | ||
136 | def test_has_wkt_property(): | |
137 | assert ( | |
138 | CRS({"init": "EPSG:4326"}) | |
139 | .to_wkt("WKT1_GDAL") | |
140 | .startswith('GEOGCS["WGS 84",DATUM') | |
141 | ) | |
142 | ||
143 | ||
144 | def test_repr(): | |
145 | assert repr(CRS({"init": "EPSG:4326"})) == ( | |
146 | "<CRS: +init=epsg:4326 +type=crs>\n" | |
147 | "Name: WGS 84\n" | |
148 | "Ellipsoid:\n" | |
149 | "- semi_major_metre: 6378137.00\n" | |
150 | "- semi_minor_metre: 6356752.31\n" | |
151 | "- inverse_flattening: 298.26\n" | |
152 | "Area of Use:\n" | |
153 | "- name: World\n" | |
154 | "- bounds: (-180.0, -90.0, 180.0, 90.0)\n" | |
155 | "Prime Meridian:\n" | |
156 | "- longitude: 0.0000\n" | |
157 | "- unit_name: degree\n" | |
158 | "- unit_conversion_factor: 0.01745329\n" | |
159 | "Axis Info:\n" | |
160 | "- Longitude[lon] (east) EPSG:9122 (degree)\n" | |
161 | "- Latitude[lat] (north) EPSG:9122 (degree)\n" | |
162 | ) | |
163 | ||
164 | ||
165 | def test_repr__long(): | |
166 | assert repr(CRS(CRS({"init": "EPSG:4326"}).to_wkt())) == ( | |
167 | '<CRS: GEOGCRS["WGS 84",DATUM["World Geodetic System 1984 ...>\n' | |
168 | "Name: WGS 84\n" | |
169 | "Ellipsoid:\n" | |
170 | "- semi_major_metre: 6378137.00\n" | |
171 | "- semi_minor_metre: 6356752.31\n" | |
172 | "- inverse_flattening: 298.26\n" | |
173 | "Area of Use:\n" | |
174 | "- name: World\n" | |
175 | "- bounds: (-180.0, -90.0, 180.0, 90.0)\n" | |
176 | "Prime Meridian:\n" | |
177 | "- longitude: 0.0000\n" | |
178 | "- unit_name: degree\n" | |
179 | "- unit_conversion_factor: 0.01745329\n" | |
180 | "Axis Info:\n" | |
181 | "- Longitude[lon] (east) EPSG:9122 (degree)\n" | |
182 | "- Latitude[lat] (north) EPSG:9122 (degree)\n" | |
183 | ) | |
184 | ||
185 | ||
186 | def test_repr__undefined(): | |
187 | assert repr( | |
188 | CRS( | |
189 | "+proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@null" | |
190 | " +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs" | |
191 | ) | |
192 | ) == ( | |
193 | "<CRS: +proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@nu ...>\n" | |
194 | "Name: unknown\n" | |
195 | "Ellipsoid:\n" | |
196 | "- semi_major_metre: 6378137.00\n" | |
197 | "- semi_minor_metre: nan\n" | |
198 | "- inverse_flattening: 0.00\n" | |
199 | "Area of Use:\n" | |
200 | "- UNDEFINED\n" | |
201 | "Prime Meridian:\n" | |
202 | "- longitude: 0.0000\n" | |
203 | "- unit_name: degree\n" | |
204 | "- unit_conversion_factor: 0.01745329\n" | |
205 | "Axis Info:\n" | |
206 | "- UNDEFINED" | |
207 | ) | |
208 | ||
209 | ||
210 | def test_dunder_str(): | |
211 | assert str(CRS({"init": "EPSG:4326"})) == CRS({"init": "EPSG:4326"}).srs | |
212 | ||
213 | ||
214 | def test_epsg(): | |
215 | assert CRS({"init": "EPSG:4326"}).to_epsg(20) == 4326 | |
216 | assert CRS({"init": "EPSG:4326"}).to_epsg() is None | |
217 | assert CRS.from_user_input(4326).to_epsg() == 4326 | |
218 | assert CRS.from_epsg(4326).to_epsg() == 4326 | |
219 | assert CRS.from_user_input("epsg:4326").to_epsg() == 4326 | |
220 | ||
221 | ||
222 | def test_datum(): | |
223 | assert CRS.from_epsg(4326).datum == CRS( | |
224 | 'DATUM["World Geodetic System 1984",' | |
225 | 'ELLIPSOID["WGS 84",6378137,298.257223563,' | |
226 | 'LENGTHUNIT["metre",1]],ID["EPSG",6326]]' | |
227 | ) | |
228 | ||
229 | ||
230 | def test_epsg__not_found(): | |
231 | assert CRS("+proj=longlat +datum=WGS84 +no_defs").to_epsg(0) is None | |
232 | assert CRS.from_string("+proj=longlat +datum=WGS84 +no_defs").to_epsg() is None | |
233 | ||
234 | ||
235 | def test_epsg__no_code_available(): | |
236 | lcc_crs = CRS.from_string( | |
237 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc " | |
238 | "+x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
239 | ) | |
240 | assert lcc_crs.to_epsg() is None | |
241 | ||
242 | ||
243 | def test_crs_OSR_equivalence(): | |
244 | crs1 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
245 | crs2 = CRS.from_string("+proj=latlong +datum=WGS84 +no_defs") | |
246 | crs3 = CRS({"init": "EPSG:4326"}) | |
247 | assert crs1 == crs2 | |
248 | # these are not equivalent in proj.4 now as one uses degrees and the othe radians | |
249 | assert crs1 == crs3 | |
250 | ||
251 | ||
252 | def test_crs_OSR_no_equivalence(): | |
253 | crs1 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
254 | crs2 = CRS.from_string("+proj=longlat +datum=NAD27 +no_defs") | |
255 | assert crs1 != crs2 | |
256 | ||
257 | ||
258 | def test_from_wkt(): | |
259 | wgs84 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
260 | from_wkt = CRS(wgs84.to_wkt()) | |
261 | assert wgs84.to_wkt() == from_wkt.to_wkt() | |
262 | ||
263 | ||
264 | def test_from_wkt_invalid(): | |
265 | with pytest.raises(CRSError): | |
266 | CRS("trash-54322") | |
267 | with pytest.raises(CRSError): | |
268 | CRS("") | |
269 | ||
270 | ||
271 | def test_from_user_input_epsg(): | |
272 | assert "+proj=longlat" in CRS.from_user_input("EPSG:4326").to_proj4(4) | |
273 | ||
274 | ||
275 | def test_from_esri_wkt(): | |
276 | projection_string = ( | |
277 | 'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",' | |
278 | 'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",' | |
279 | 'SPHEROID["GRS_1980",6378137.0,298.257222101]],' | |
280 | 'PRIMEM["Greenwich",0.0],' | |
281 | 'UNIT["Degree",0.0174532925199433]],' | |
282 | 'PROJECTION["Albers"],' | |
283 | 'PARAMETER["false_easting",0.0],' | |
284 | 'PARAMETER["false_northing",0.0],' | |
285 | 'PARAMETER["central_meridian",-96.0],' | |
286 | 'PARAMETER["standard_parallel_1",29.5],' | |
287 | 'PARAMETER["standard_parallel_2",45.5],' | |
288 | 'PARAMETER["latitude_of_origin",23.0],' | |
289 | 'UNIT["Meter",1.0],' | |
290 | 'VERTCS["NAVD_1988",' | |
291 | 'VDATUM["North_American_Vertical_Datum_1988"],' | |
292 | 'PARAMETER["Vertical_Shift",0.0],' | |
293 | 'PARAMETER["Direction",1.0],UNIT["Centimeter",0.01]]]' | |
294 | ) | |
295 | proj_crs_str = CRS.from_string(projection_string) | |
296 | proj_crs_wkt = CRS(projection_string) | |
297 | assert proj_crs_str.to_proj4() == proj_crs_wkt.to_proj4() | |
298 | assert proj_crs_str.to_proj4(4) == ( | |
299 | "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 " | |
300 | "+lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs" | |
301 | ) | |
302 | ||
303 | ||
304 | def test_compound_crs(): | |
305 | wkt = """COMPD_CS["unknown",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],VERT_CS["unknown",VERT_DATUM["unknown",2005],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Up",UP]]]""" | |
306 | assert CRS(wkt).to_wkt("WKT1_GDAL").startswith('COMPD_CS["unknown",GEOGCS["WGS 84"') | |
307 | ||
308 | ||
309 | def test_ellipsoid(): | |
310 | crs1 = CRS.from_epsg(4326) | |
311 | assert "{:.3f}".format(crs1.ellipsoid.inverse_flattening) == "298.257" | |
312 | assert "{:.3f}".format(crs1.ellipsoid.semi_major_metre) == "6378137.000" | |
313 | assert "{:.3f}".format(crs1.ellipsoid.semi_minor_metre) == "6356752.314" | |
314 | ||
315 | ||
316 | def test_area_of_use(): | |
317 | crs1 = CRS.from_epsg(4326) | |
318 | assert crs1.area_of_use.bounds == (-180.0, -90.0, 180.0, 90.0) | |
319 | assert crs1.area_of_use.name == "World" | |
320 | ||
321 | ||
322 | def test_from_user_input_custom_crs_class(): | |
323 | class CustomCRS(object): | |
324 | def to_wkt(self): | |
325 | return CRS.from_epsg(4326).to_wkt() | |
326 | ||
327 | assert CRS.from_user_input(CustomCRS()) == CRS.from_epsg(4326) |
0 | import os | |
1 | import shutil | |
2 | import tempfile | |
3 | import unittest | |
4 | from contextlib import contextmanager | |
5 | ||
6 | import pytest | |
7 | from mock import patch | |
8 | from pyproj.datadir import DataDirError, append_data_dir, get_data_dir, set_data_dir | |
9 | ||
10 | ||
11 | def create_projdb(tmpdir): | |
12 | with open(os.path.join(tmpdir, "proj.db"), "w") as pjdb: | |
13 | pjdb.write("DUMMY proj.db") | |
14 | ||
15 | ||
16 | @contextmanager | |
17 | def proj_env(): | |
18 | """ | |
19 | Ensure environment variable the same at the end of the test. | |
20 | """ | |
21 | proj_lib = os.environ.get("PROJ_LIB") | |
22 | try: | |
23 | yield | |
24 | finally: | |
25 | if proj_lib: | |
26 | os.environ["PROJ_LIB"] = proj_lib | |
27 | ||
28 | ||
29 | @contextmanager | |
30 | def temporary_directory(): | |
31 | """ | |
32 | Get a temporary directory | |
33 | """ | |
34 | temp_dir = tempfile.mkdtemp() | |
35 | try: | |
36 | yield temp_dir | |
37 | finally: | |
38 | shutil.rmtree(temp_dir) | |
39 | ||
40 | ||
41 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
42 | def test_get_data_dir__missing(): | |
43 | with proj_env(), pytest.raises(DataDirError), patch( | |
44 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
45 | ), patch("pyproj.datadir.find_executable", return_value=None): | |
46 | set_data_dir(None) | |
47 | os.environ.pop("PROJ_LIB", None) | |
48 | get_data_dir() | |
49 | ||
50 | ||
51 | def test_get_data_dir__from_user(): | |
52 | with proj_env(), temporary_directory() as tmpdir, temporary_directory() as tmpdir_env: | |
53 | create_projdb(tmpdir) | |
54 | os.environ["PROJ_LIB"] = tmpdir_env | |
55 | create_projdb(tmpdir_env) | |
56 | set_data_dir(tmpdir) | |
57 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
58 | os.makedirs(internal_proj_dir) | |
59 | create_projdb(internal_proj_dir) | |
60 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
61 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
62 | assert get_data_dir() == tmpdir | |
63 | ||
64 | ||
65 | def test_get_data_dir__internal(): | |
66 | with proj_env(), temporary_directory() as tmpdir: | |
67 | set_data_dir(None) | |
68 | os.environ["PROJ_LIB"] = tmpdir | |
69 | create_projdb(tmpdir) | |
70 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
71 | os.makedirs(internal_proj_dir) | |
72 | create_projdb(internal_proj_dir) | |
73 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
74 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
75 | assert get_data_dir() == internal_proj_dir | |
76 | ||
77 | ||
78 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
79 | def test_get_data_dir__from_env_var(): | |
80 | with proj_env(), temporary_directory() as tmpdir, patch( | |
81 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
82 | ): | |
83 | set_data_dir(None) | |
84 | os.environ["PROJ_LIB"] = tmpdir | |
85 | create_projdb(tmpdir) | |
86 | assert get_data_dir() == tmpdir | |
87 | ||
88 | ||
89 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
90 | def test_get_data_dir__from_env_var__multiple(): | |
91 | with proj_env(), temporary_directory() as tmpdir, patch( | |
92 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
93 | ): | |
94 | set_data_dir(None) | |
95 | os.environ["PROJ_LIB"] = os.pathsep.join([tmpdir, tmpdir, tmpdir]) | |
96 | create_projdb(tmpdir) | |
97 | assert get_data_dir() == os.pathsep.join([tmpdir, tmpdir, tmpdir]) | |
98 | ||
99 | ||
100 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
101 | def test_get_data_dir__from_path(): | |
102 | with proj_env(), temporary_directory() as tmpdir, patch( | |
103 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
104 | ), patch("pyproj.datadir.find_executable") as find_exe: | |
105 | set_data_dir(None) | |
106 | os.environ.pop("PROJ_LIB", None) | |
107 | find_exe.return_value = os.path.join(tmpdir, "bin", "proj") | |
108 | proj_dir = os.path.join(tmpdir, "share", "proj") | |
109 | os.makedirs(proj_dir) | |
110 | create_projdb(proj_dir) | |
111 | assert get_data_dir() == proj_dir | |
112 | ||
113 | ||
114 | def test_append_data_dir__internal(): | |
115 | with proj_env(), temporary_directory() as tmpdir: | |
116 | set_data_dir(None) | |
117 | os.environ["PROJ_LIB"] = tmpdir | |
118 | create_projdb(tmpdir) | |
119 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
120 | os.makedirs(internal_proj_dir) | |
121 | create_projdb(internal_proj_dir) | |
122 | extra_datadir = str(os.path.join(tmpdir, "extra_datumgrids")) | |
123 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
124 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
125 | append_data_dir(extra_datadir) | |
126 | assert get_data_dir() == os.pathsep.join([internal_proj_dir, extra_datadir]) |
0 | import subprocess, sys | |
0 | from numpy.testing import assert_almost_equal | |
1 | ||
1 | 2 | from pyproj import Proj, transform |
2 | 3 | |
3 | p1 = Proj(proj="latlong", datum="WGS84") | |
4 | s_1 = -111.5 | |
5 | s_2 = 45.25919444444 | |
6 | p2 = Proj(proj="utm", zone=10, datum="NAD27") | |
7 | cstr = ( | |
8 | 'echo "-111.5 45.25919444444" | cs2cs +proj=latlong +datum=WGS84 +to ' | |
9 | + '+proj=utm +zone=10 +datum=NAD27 -f "%.12f"' | |
10 | ) | |
11 | sys.stdout.write("test datum shift (requires cs2cs proj4 command line tool)\n") | |
12 | p = subprocess.Popen(cstr, shell=True, stdout=subprocess.PIPE, stderr=subprocess.STDOUT) | |
13 | pout = p.stdout.readlines()[0].decode("ascii") | |
14 | x1, y1, z1 = tuple([float(x) for x in pout.split()]) | |
15 | x2, y2 = transform(p1, p2, s_1, s_2) | |
16 | sys.stdout.write("is (%s,%s) = (%s,%s)?\n" % (x1, y1, x2, y2)) | |
17 | print("%12.3f %12.3f" % (x1, y1) == "%12.3f %12.3f" % (x2, y2)) | |
4 | ||
5 | def test_datum(): | |
6 | p1 = Proj(proj="latlong", datum="WGS84") | |
7 | s_1 = -111.5 | |
8 | s_2 = 45.25919444444 | |
9 | p2 = Proj(proj="utm", zone=10, datum="NAD27") | |
10 | x2, y2 = transform(p1, p2, s_1, s_2) | |
11 | assert_almost_equal((x2, y2), (1402291.0833290431, 5076289.591846835)) |
0 | import pytest | |
1 | from numpy.testing import assert_almost_equal | |
2 | ||
3 | from pyproj import Proj, transform | |
4 | ||
5 | # illustrates the use of the transform function to | |
6 | # perform coordinate transformations with datum shifts. | |
7 | # | |
8 | # This example is from Roberto Vidmar | |
9 | # | |
10 | # Test point is Trieste, Molo Sartorio | |
11 | # | |
12 | # This data come from the Istituto Geografico Militare (IGM), as well as | |
13 | # the 7 parameters to transform from Gauss-Boaga (our reference frame) | |
14 | # to WGS84 | |
15 | # | |
16 | # WGS84 Lat: 45d38'49.879" (45.647188611) | |
17 | # WGS84 Lon: 13d45'34.397" (13.759554722) | |
18 | # WGS84 z: 52.80; | |
19 | # UTM 33: 403340.97 5055597.17 | |
20 | # GB: 2423346.99 5055619.87 | |
21 | ||
22 | UTM_x = 403340.9672367854 | |
23 | UTM_y = 5055597.175553089 | |
24 | GB_x = 2423346.99 | |
25 | GB_y = 5055619.87 | |
26 | WGS84_lat = 45.647188611 # Degrees | |
27 | WGS84_lon = 13.759554722 # Degrees | |
28 | UTM_z = WGS84_z = 52.8 # Ellipsoidical height in meters | |
29 | WGS84_PROJ = Proj(proj="latlong", datum="WGS84") | |
30 | UTM_33_PROJ = Proj(proj="utm", zone="33") | |
31 | GAUSSSB_PROJ = Proj( | |
32 | init="epsg:3004", towgs84="-122.74,-34.27,-22.83,-1.884,-3.400,-3.030,-15.62" | |
33 | ) | |
34 | print("proj4 library version = ", WGS84_PROJ.proj_version) | |
35 | ||
36 | ||
37 | def test_shift_wgs84_to_utm33(): | |
38 | xutm33, yutm33, zutm33 = transform( | |
39 | WGS84_PROJ, UTM_33_PROJ, WGS84_lon, WGS84_lat, WGS84_z | |
40 | ) | |
41 | assert_almost_equal((xutm33, yutm33, zutm33), (UTM_x, UTM_y, UTM_z)) | |
42 | ||
43 | ||
44 | def test_shift_utm33_to_wgs84(): | |
45 | back_lon, back_lat, back_z = transform(UTM_33_PROJ, WGS84_PROJ, UTM_x, UTM_y, UTM_z) | |
46 | assert_almost_equal((back_lon, back_lat, back_z), (WGS84_lon, WGS84_lat, WGS84_z)) | |
47 | ||
48 | ||
49 | def test_shift_wgs84_to_gaussb_no_ellisoidal_height(): | |
50 | xgb, ygb, zgb = transform(WGS84_PROJ, GAUSSSB_PROJ, WGS84_lon, WGS84_lat, 0) | |
51 | assert_almost_equal((xgb, ygb, zgb), (GB_x, 5055619.899, 0), decimal=2) | |
52 | ||
53 | ||
54 | def test_shift_gaussb_to_wgs84_no_ellisoidal_height(): | |
55 | back_lon, back_lat, back_z = transform(GAUSSSB_PROJ, WGS84_PROJ, GB_x, GB_y, 0) | |
56 | assert_almost_equal( | |
57 | (back_lon, back_lat, back_z), (WGS84_lon, WGS84_lat, 0), decimal=3 | |
58 | ) |
0 | """ | |
1 | This is a wrapper for the doctests in lib/pyproj/__init__.py so that | |
2 | pytest can conveniently run all the tests in a single command line. | |
3 | """ | |
4 | ||
5 | import pyproj | |
6 | ||
7 | ||
8 | def test_doctests(): | |
9 | failure_count = pyproj.test() | |
10 | # if the below line fails, doctests have failed | |
11 | assert ( | |
12 | failure_count == 0 | |
13 | ), "{0} of the doctests in " "lib/pyproj/__init__.py failed".format(failure_count) |
0 | import pytest | |
1 | ||
2 | from pyproj import CRS, Proj | |
3 | from pyproj.exceptions import CRSError, ProjError | |
4 | ||
5 | ||
6 | def test_proj_exception(): | |
7 | with pytest.raises(ProjError, match="Internal Proj Error"): | |
8 | Proj("+proj=bobbyjoe") | |
9 | ||
10 | ||
11 | def test_crs_exception(): | |
12 | with pytest.raises(CRSError, match="Internal Proj Error"): | |
13 | CRS("+proj=bobbyjoe") |
0 | import os | |
1 | import pickle | |
2 | import shutil | |
3 | import tempfile | |
4 | from contextlib import contextmanager | |
5 | ||
6 | from numpy.testing import assert_almost_equal | |
7 | ||
8 | from pyproj import Geod | |
9 | ||
10 | ||
11 | @contextmanager | |
12 | def temporary_directory(): | |
13 | """ | |
14 | Get a temporary directory | |
15 | """ | |
16 | temp_dir = tempfile.mkdtemp() | |
17 | try: | |
18 | yield temp_dir | |
19 | finally: | |
20 | shutil.rmtree(temp_dir) | |
21 | ||
22 | ||
23 | def test_geod_inverse_transform(): | |
24 | gg = Geod(ellps="clrk66") | |
25 | lat1pt = 42.0 + (15.0 / 60.0) | |
26 | lon1pt = -71.0 - (7.0 / 60.0) | |
27 | lat2pt = 45.0 + (31.0 / 60.0) | |
28 | lon2pt = -123.0 - (41.0 / 60.0) | |
29 | """ | |
30 | distance between boston and portland, clrk66: | |
31 | -66.531 75.654 4164192.708 | |
32 | distance between boston and portland, WGS84: | |
33 | -66.530 75.654 4164074.239 | |
34 | testing pickling of Geod instance | |
35 | distance between boston and portland, clrk66 (from pickle): | |
36 | -66.531 75.654 4164192.708 | |
37 | distance between boston and portland, WGS84 (from pickle): | |
38 | -66.530 75.654 4164074.239 | |
39 | inverse transform | |
40 | from proj.4 invgeod: | |
41 | b'-66.531\t75.654\t4164192.708\n' | |
42 | ||
43 | """ | |
44 | print("from pyproj.Geod.inv:") | |
45 | az12, az21, dist = gg.inv(lon1pt, lat1pt, lon2pt, lat2pt) | |
46 | assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) | |
47 | ||
48 | print("forward transform") | |
49 | print("from proj.4 geod:") | |
50 | endlon, endlat, backaz = gg.fwd(lon1pt, lat1pt, az12, dist) | |
51 | assert_almost_equal((endlon, endlat, backaz), (-123.683, 45.517, 75.654), decimal=3) | |
52 | print("intermediate points:") | |
53 | print("from geod with +lat_1,+lon_1,+lat_2,+lon_2,+n_S:") | |
54 | npts = 4 | |
55 | lonlats = gg.npts(lon1pt, lat1pt, lon2pt, lat2pt, npts) | |
56 | lonprev = lon1pt | |
57 | latprev = lat1pt | |
58 | print(dist / (npts + 1)) | |
59 | print("%6.3f %7.3f" % (lat1pt, lon1pt)) | |
60 | result_dists = ( | |
61 | (-66.53059478766238, 106.79071710136431, 832838.5416198927), | |
62 | (-73.20928289863558, 99.32289055927389, 832838.5416198935), | |
63 | (-80.67710944072617, 91.36325611787134, 832838.5416198947), | |
64 | (-88.63674388212858, 83.32809401477382, 832838.5416198922), | |
65 | ) | |
66 | for (lon, lat), (res12, res21, resdist) in zip(lonlats, result_dists): | |
67 | az12, az21, dist = gg.inv(lonprev, latprev, lon, lat) | |
68 | assert_almost_equal((az12, az21, dist), (res12, res21, resdist)) | |
69 | latprev = lat | |
70 | lonprev = lon | |
71 | az12, az21, dist = gg.inv(lonprev, latprev, lon2pt, lat2pt) | |
72 | assert_almost_equal( | |
73 | (lat2pt, lon2pt, dist), (45.517, -123.683, 832838.542), decimal=3 | |
74 | ) | |
75 | ||
76 | ||
77 | def test_geod_cities(): | |
78 | # specify the lat/lons of some cities. | |
79 | boston_lat = 42.0 + (15.0 / 60.0) | |
80 | boston_lon = -71.0 - (7.0 / 60.0) | |
81 | portland_lat = 45.0 + (31.0 / 60.0) | |
82 | portland_lon = -123.0 - (41.0 / 60.0) | |
83 | g1 = Geod(ellps="clrk66") | |
84 | g2 = Geod(ellps="WGS84") | |
85 | az12, az21, dist = g1.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
86 | print("distance between boston and portland, clrk66:") | |
87 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
88 | assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) | |
89 | print("distance between boston and portland, WGS84:") | |
90 | az12, az21, dist = g2.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
91 | assert_almost_equal((az12, az21, dist), (-66.530, 75.654, 4164074.239), decimal=3) | |
92 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
93 | print("testing pickling of Geod instance") | |
94 | with temporary_directory() as tmpdir: | |
95 | with open(os.path.join(tmpdir, "geod1.pickle"), "wb") as gp1w: | |
96 | pickle.dump(g1, gp1w, -1) | |
97 | with open(os.path.join(tmpdir, "geod2.pickle"), "wb") as gp2w: | |
98 | pickle.dump(g2, gp2w, -1) | |
99 | with open(os.path.join(tmpdir, "geod1.pickle"), "rb") as gp1: | |
100 | g3 = pickle.load(gp1) | |
101 | with open(os.path.join(tmpdir, "geod2.pickle"), "rb") as gp2: | |
102 | g4 = pickle.load(gp2) | |
103 | az12, az21, dist = g3.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
104 | assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) | |
105 | print("distance between boston and portland, clrk66 (from pickle):") | |
106 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
107 | az12, az21, dist = g4.inv(boston_lon, boston_lat, portland_lon, portland_lat) | |
108 | print("distance between boston and portland, WGS84 (from pickle):") | |
109 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
110 | assert_almost_equal((az12, az21, dist), (-66.530, 75.654, 4164074.239), decimal=3) | |
111 | g3 = Geod("+ellps=clrk66") # proj4 style init string | |
112 | print("inverse transform") | |
113 | lat1pt = 42.0 + (15.0 / 60.0) | |
114 | lon1pt = -71.0 - (7.0 / 60.0) | |
115 | lat2pt = 45.0 + (31.0 / 60.0) | |
116 | lon2pt = -123.0 - (41.0 / 60.0) | |
117 | az12, az21, dist = g3.inv(lon1pt, lat1pt, lon2pt, lat2pt) | |
118 | print("%7.3f %6.3f %12.3f" % (az12, az21, dist)) | |
119 | assert_almost_equal((az12, az21, dist), (-66.531, 75.654, 4164192.708), decimal=3) |
0 | """run test.py first!""" | |
1 | import os | |
2 | import pickle | |
3 | import shutil | |
4 | import tempfile | |
5 | import time | |
6 | from contextlib import contextmanager | |
7 | ||
8 | import numpy | |
9 | from numpy.testing import assert_allclose | |
10 | ||
11 | from pyproj import Proj | |
12 | ||
13 | ||
14 | @contextmanager | |
15 | def temporary_directory(): | |
16 | """ | |
17 | Get a temporary directory | |
18 | """ | |
19 | temp_dir = tempfile.mkdtemp() | |
20 | try: | |
21 | yield temp_dir | |
22 | finally: | |
23 | shutil.rmtree(temp_dir) | |
24 | ||
25 | ||
26 | def test_pickle(): | |
27 | nx = 349 | |
28 | ny = 277 | |
29 | dx = 32463.41 | |
30 | dy = dx | |
31 | print("do it again, from pickled instance ...") | |
32 | # find 4 lon/lat crnrs of AWIPS grid 221. | |
33 | llcrnrx = 0.0 | |
34 | llcrnry = 0.0 | |
35 | urcrnrx = dx * (nx - 1) | |
36 | urcrnry = dy * (ny - 1) | |
37 | dx = (urcrnrx - llcrnrx) / (nx - 1) | |
38 | dy = (urcrnry - llcrnry) / (ny - 1) | |
39 | x = llcrnrx + dx * numpy.indices((ny, nx), "f")[1, :, :] | |
40 | y = llcrnry + dy * numpy.indices((ny, nx), "f")[0, :, :] | |
41 | awips221_pre_pickle = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
42 | with temporary_directory() as tmpdir: | |
43 | with open(os.path.join(tmpdir, "test.pickle"), "wb") as pfh: | |
44 | pickle.dump(awips221_pre_pickle, pfh, -1) | |
45 | with open(os.path.join(tmpdir, "test.pickle"), "rb") as prh: | |
46 | awips221 = pickle.load(prh) | |
47 | t1 = time.clock() | |
48 | lons, lats = awips221(x, y, inverse=True) | |
49 | t2 = time.clock() | |
50 | print("compute lats/lons for all points on AWIPS 221 grid (%sx%s)" % (nx, ny)) | |
51 | print("max/min lons in radians") | |
52 | print( | |
53 | numpy.minimum.reduce(numpy.ravel(lons)), numpy.maximum.reduce(numpy.ravel(lons)) | |
54 | ) | |
55 | print("max/min lats in radians") | |
56 | print( | |
57 | numpy.minimum.reduce(numpy.ravel(lats)), numpy.maximum.reduce(numpy.ravel(lats)) | |
58 | ) | |
59 | print("took", t2 - t1, "secs") | |
60 | # reverse transformation. | |
61 | t1 = time.clock() | |
62 | xx, yy = awips221(lons, lats) | |
63 | t2 = time.clock() | |
64 | print("max abs error for x") | |
65 | max_abs_err_x = numpy.maximum.reduce(numpy.fabs(numpy.ravel(x - xx))) | |
66 | print(max_abs_err_x) | |
67 | assert_allclose(max_abs_err_x, 0, atol=1e-4) | |
68 | print("max abs error for y") | |
69 | max_abs_err_y = numpy.maximum.reduce(numpy.fabs(numpy.ravel(y - yy))) | |
70 | print(max_abs_err_y) | |
71 | assert_allclose(max_abs_err_y, 0, atol=1e-4) | |
72 | print("took", t2 - t1, "secs") |
0 | import numpy | |
1 | from numpy.testing import assert_allclose | |
2 | ||
0 | 3 | from pyproj import Proj, transform |
1 | import numpy as N | |
2 | 4 | |
3 | # convert awips221 grid to awips218 coordinate system | |
4 | # (grids defined at http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html) | |
5 | nx = 614 | |
6 | ny = 428 | |
7 | dx = 12190.58 | |
8 | dy = dx | |
9 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
10 | print("proj4 library version = ", awips221.proj_version) | |
11 | llcrnrx, llcrnry = awips221(-145.5, 1) | |
12 | awips221 = Proj( | |
13 | proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107, x_0=-llcrnrx, y_0=-llcrnry | |
14 | ) | |
15 | print(awips221(-145.5, 1), "(should be close to zero)") | |
16 | awips218 = Proj(proj="lcc", R=6371200, lat_1=25, lat_2=25, lon_0=-95) | |
17 | llcrnrx, llcrnry = awips218(-133.459, 12.19) | |
18 | awips218 = Proj( | |
19 | proj="lcc", R=6371200, lat_1=25, lat_2=25, lon_0=-95, x_0=-llcrnrx, y_0=-llcrnry | |
20 | ) | |
21 | print(awips218(-133.459, 12.19), "(should close to zero)") | |
22 | x1 = dx * N.indices((ny, nx), "f")[1, :, :] | |
23 | y1 = dy * N.indices((ny, nx), "f")[0, :, :] | |
24 | print("max/min x and y for awips218 grid") | |
25 | print(N.minimum.reduce(N.ravel(x1)), N.maximum.reduce(N.ravel(x1))) | |
26 | print(N.minimum.reduce(N.ravel(y1)), N.maximum.reduce(N.ravel(y1))) | |
27 | x2, y2 = transform(awips218, awips221, x1, y1) | |
28 | print("max/min x and y for awips218 grid in awips221 coordinates") | |
29 | print(N.minimum.reduce(N.ravel(x2)), N.maximum.reduce(N.ravel(x2))) | |
30 | print(N.minimum.reduce(N.ravel(y2)), N.maximum.reduce(N.ravel(y2))) | |
31 | x3, y3 = transform(awips221, awips218, x2, y2) | |
32 | print("error for reverse transformation back to awips218 coords") | |
33 | print("(should be close to zero)") | |
34 | print(N.minimum.reduce(N.ravel(x3 - x1)), N.maximum.reduce(N.ravel(x3 - x1))) | |
35 | print(N.minimum.reduce(N.ravel(y3 - y1)), N.maximum.reduce(N.ravel(y3 - y1))) | |
5 | ||
6 | def test_transform(): | |
7 | # convert awips221 grid to awips218 coordinate system | |
8 | # (grids defined at http://www.nco.ncep.noaa.gov/pmb/docs/on388/tableb.html) | |
9 | nx = 614 | |
10 | ny = 428 | |
11 | dx = 12190.58 | |
12 | dy = dx | |
13 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
14 | print("proj4 library version = ", awips221.proj_version) | |
15 | llcrnrx, llcrnry = awips221(-145.5, 1) | |
16 | awips221 = Proj( | |
17 | proj="lcc", | |
18 | R=6371200, | |
19 | lat_1=50, | |
20 | lat_2=50, | |
21 | lon_0=-107, | |
22 | x_0=-llcrnrx, | |
23 | y_0=-llcrnry, | |
24 | ) | |
25 | assert_allclose(awips221(-145.5, 1), (0, 0), atol=1e-4) | |
26 | awips218 = Proj(proj="lcc", R=6371200, lat_1=25, lat_2=25, lon_0=-95) | |
27 | llcrnrx, llcrnry = awips218(-133.459, 12.19) | |
28 | awips218 = Proj( | |
29 | proj="lcc", R=6371200, lat_1=25, lat_2=25, lon_0=-95, x_0=-llcrnrx, y_0=-llcrnry | |
30 | ) | |
31 | assert_allclose(awips218(-133.459, 12.19), (0, 0), atol=1e-4) | |
32 | x1 = dx * numpy.indices((ny, nx), "f")[1, :, :] | |
33 | y1 = dy * numpy.indices((ny, nx), "f")[0, :, :] | |
34 | print("max/min x and y for awips218 grid") | |
35 | print(numpy.minimum.reduce(numpy.ravel(x1)), numpy.maximum.reduce(numpy.ravel(x1))) | |
36 | print(numpy.minimum.reduce(numpy.ravel(y1)), numpy.maximum.reduce(numpy.ravel(y1))) | |
37 | x2, y2 = transform(awips218, awips221, x1, y1) | |
38 | print("max/min x and y for awips218 grid in awips221 coordinates") | |
39 | print(numpy.minimum.reduce(numpy.ravel(x2)), numpy.maximum.reduce(numpy.ravel(x2))) | |
40 | print(numpy.minimum.reduce(numpy.ravel(y2)), numpy.maximum.reduce(numpy.ravel(y2))) | |
41 | x3, y3 = transform(awips221, awips218, x2, y2) | |
42 | print("error for reverse transformation back to awips218 coords") | |
43 | print("(should be close to zero)") | |
44 | assert_allclose(numpy.minimum.reduce(numpy.ravel(x3 - x1)), 0, atol=1e-4) | |
45 | assert_allclose(numpy.maximum.reduce(numpy.ravel(x3 - x1)), 0, atol=1e-4) | |
46 | assert_allclose(numpy.minimum.reduce(numpy.ravel(y3 - y1)), 0, atol=1e-4) | |
47 | assert_allclose(numpy.maximum.reduce(numpy.ravel(y3 - y1)), 0, atol=1e-4) |
0 | import numpy as np | |
1 | import pytest | |
2 | from numpy.testing import assert_almost_equal, assert_equal | |
3 | ||
4 | import pyproj | |
5 | from pyproj import Proj, Transformer, itransform, transform | |
6 | from pyproj.exceptions import ProjError | |
7 | ||
8 | ||
9 | def test_tranform_wgs84_to_custom(): | |
10 | custom_proj = pyproj.Proj( | |
11 | "+proj=geos +lon_0=0.000000 +lat_0=0 +h=35807.414063" | |
12 | " +a=6378.169000 +b=6356.583984" | |
13 | ) | |
14 | wgs84 = pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
15 | lat, lon = 51.04715, 3.23406 | |
16 | xx, yy = pyproj.transform(wgs84, custom_proj, lon, lat) | |
17 | assert "{:.3f} {:.3f}".format(xx, yy) == "212.623 4604.975" | |
18 | ||
19 | ||
20 | def test_transform_wgs84_to_alaska(): | |
21 | lat_lon_proj = pyproj.Proj(init="epsg:4326", preserve_units=False) | |
22 | alaska_aea_proj = pyproj.Proj(init="epsg:2964", preserve_units=False) | |
23 | test = (-179.72638, 49.752533) | |
24 | xx, yy = pyproj.transform(lat_lon_proj, alaska_aea_proj, *test) | |
25 | assert "{:.3f} {:.3f}".format(xx, yy) == "-1824924.495 330822.800" | |
26 | ||
27 | ||
28 | def test_illegal_transformation(): | |
29 | # issue 202 | |
30 | p1 = pyproj.Proj(init="epsg:4326") | |
31 | p2 = pyproj.Proj(init="epsg:3857") | |
32 | xx, yy = pyproj.transform( | |
33 | p1, p2, (-180, -180, 180, 180, -180), (-90, 90, 90, -90, -90) | |
34 | ) | |
35 | assert np.all(np.isinf(xx)) | |
36 | assert np.all(np.isinf(yy)) | |
37 | try: | |
38 | xx, yy = pyproj.transform( | |
39 | p1, p2, (-180, -180, 180, 180, -180), (-90, 90, 90, -90, -90), errcheck=True | |
40 | ) | |
41 | assert_equal(None, "Should throw an exception when errcheck=True") | |
42 | except ProjError: | |
43 | pass | |
44 | ||
45 | ||
46 | def test_lambert_conformal_transform(): | |
47 | # issue 207 | |
48 | Midelt = pyproj.Proj(init="epsg:26191") | |
49 | WGS84 = pyproj.Proj(init="epsg:4326") | |
50 | ||
51 | E = 567623.931 | |
52 | N = 256422.787 | |
53 | h = 1341.467 | |
54 | ||
55 | Long1, Lat1, H1 = pyproj.transform(Midelt, WGS84, E, N, h, radians=False) | |
56 | assert_almost_equal((Long1, Lat1, H1), (-4.6753456, 32.902199, 1341.467), decimal=5) | |
57 | ||
58 | ||
59 | def test_equivalent_crs(): | |
60 | transformer = Transformer.from_crs("epsg:4326", 4326, skip_equivalent=True) | |
61 | assert transformer._transformer.projections_equivalent | |
62 | assert transformer._transformer.projections_exact_same | |
63 | assert transformer._transformer.skip_equivalent | |
64 | ||
65 | ||
66 | def test_equivalent_crs__disabled(): | |
67 | transformer = Transformer.from_crs("epsg:4326", 4326) | |
68 | assert not transformer._transformer.skip_equivalent | |
69 | assert transformer._transformer.projections_equivalent | |
70 | assert transformer._transformer.projections_exact_same | |
71 | ||
72 | ||
73 | def test_equivalent_crs__different(): | |
74 | transformer = Transformer.from_crs("epsg:4326", 3857, skip_equivalent=True) | |
75 | assert transformer._transformer.skip_equivalent | |
76 | assert not transformer._transformer.projections_equivalent | |
77 | assert not transformer._transformer.projections_exact_same | |
78 | ||
79 | ||
80 | def test_equivalent_proj(): | |
81 | transformer = Transformer.from_proj( | |
82 | "+init=epsg:4326", pyproj.Proj(4326).crs.to_proj4(), skip_equivalent=True | |
83 | ) | |
84 | assert transformer._transformer.skip_equivalent | |
85 | assert transformer._transformer.projections_equivalent | |
86 | assert not transformer._transformer.projections_exact_same | |
87 | ||
88 | ||
89 | def test_equivalent_proj__disabled(): | |
90 | transformer = Transformer.from_proj(3857, pyproj.Proj(3857).crs.to_proj4()) | |
91 | assert not transformer._transformer.skip_equivalent | |
92 | assert not transformer._transformer.projections_equivalent | |
93 | assert not transformer._transformer.projections_exact_same | |
94 | ||
95 | ||
96 | def test_equivalent_proj__different(): | |
97 | transformer = Transformer.from_proj(3857, 4326, skip_equivalent=True) | |
98 | assert transformer._transformer.skip_equivalent | |
99 | assert not transformer._transformer.projections_equivalent | |
100 | assert not transformer._transformer.projections_exact_same | |
101 | ||
102 | ||
103 | def test_equivalent_pipeline(): | |
104 | transformer = Transformer.from_pipeline( | |
105 | "+proj=pipeline +step +proj=longlat +ellps=WGS84 +step " | |
106 | "+proj=unitconvert +xy_in=rad +xy_out=deg" | |
107 | ) | |
108 | assert not transformer._transformer.skip_equivalent | |
109 | assert not transformer._transformer.projections_equivalent | |
110 | assert not transformer._transformer.projections_exact_same | |
111 | ||
112 | ||
113 | def test_4d_transform(): | |
114 | transformer = Transformer.from_pipeline("+init=ITRF2008:ITRF2000") | |
115 | assert_almost_equal( | |
116 | transformer.transform( | |
117 | xx=3513638.19380, yy=778956.45250, zz=5248216.46900, tt=2008.75 | |
118 | ), | |
119 | (3513638.1999428216, 778956.4532640711, 5248216.453456361, 2008.75), | |
120 | ) | |
121 | ||
122 | ||
123 | def test_2d_with_time_transform(): | |
124 | transformer = Transformer.from_pipeline("+init=ITRF2008:ITRF2000") | |
125 | assert_almost_equal( | |
126 | transformer.transform(xx=3513638.19380, yy=778956.45250, tt=2008.75), | |
127 | (3513638.1999428216, 778956.4532640711, 2008.75), | |
128 | ) | |
129 | ||
130 | ||
131 | def test_4d_transform_crs_obs1(): | |
132 | transformer = Transformer.from_proj(7789, 8401) | |
133 | assert_almost_equal( | |
134 | transformer.transform( | |
135 | xx=3496737.2679, yy=743254.4507, zz=5264462.9620, tt=2019.0 | |
136 | ), | |
137 | (3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0), | |
138 | ) | |
139 | ||
140 | ||
141 | def test_4d_transform_orginal_crs_obs1(): | |
142 | assert_almost_equal( | |
143 | transform(7789, 8401, x=3496737.2679, y=743254.4507, z=5264462.9620, tt=2019.0), | |
144 | (3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0), | |
145 | ) | |
146 | ||
147 | ||
148 | def test_4d_transform_crs_obs2(): | |
149 | transformer = Transformer.from_proj(4896, 7930) | |
150 | assert_almost_equal( | |
151 | transformer.transform( | |
152 | xx=3496737.2679, yy=743254.4507, zz=5264462.9620, tt=2019.0 | |
153 | ), | |
154 | (3496737.7857162016, 743254.0394113371, 5264462.643659916, 2019.0), | |
155 | ) | |
156 | ||
157 | ||
158 | def test_2d_with_time_transform_crs_obs2(): | |
159 | transformer = Transformer.from_proj(4896, 7930) | |
160 | assert_almost_equal( | |
161 | transformer.transform(xx=3496737.2679, yy=743254.4507, tt=2019.0), | |
162 | (3496737.4105305015, 743254.1014318303, 2019.0), | |
163 | ) | |
164 | ||
165 | ||
166 | def test_2d_with_time_transform_original_crs_obs2(): | |
167 | assert_almost_equal( | |
168 | transform(4896, 7930, x=3496737.2679, y=743254.4507, tt=2019.0), | |
169 | (3496737.4105305015, 743254.1014318303, 2019.0), | |
170 | ) | |
171 | ||
172 | ||
173 | def test_4d_itransform(): | |
174 | transformer = Transformer.from_pipeline("+init=ITRF2008:ITRF2000") | |
175 | assert_almost_equal( | |
176 | list( | |
177 | transformer.itransform( | |
178 | [(3513638.19380, 778956.45250, 5248216.46900, 2008.75)] | |
179 | ) | |
180 | ), | |
181 | [(3513638.1999428216, 778956.4532640711, 5248216.453456361, 2008.75)], | |
182 | ) | |
183 | ||
184 | ||
185 | def test_3d_time_itransform(): | |
186 | transformer = Transformer.from_pipeline("+init=ITRF2008:ITRF2000") | |
187 | assert_almost_equal( | |
188 | list( | |
189 | transformer.itransform( | |
190 | [(3513638.19380, 778956.45250, 2008.75)], time_3rd=True | |
191 | ) | |
192 | ), | |
193 | [(3513638.1999428216, 778956.4532640711, 2008.75)], | |
194 | ) | |
195 | ||
196 | ||
197 | def test_4d_itransform_orginal_crs_obs1(): | |
198 | assert_almost_equal( | |
199 | list( | |
200 | itransform(7789, 8401, [(3496737.2679, 743254.4507, 5264462.9620, 2019.0)]) | |
201 | ), | |
202 | [(3496737.757717311, 743253.9940103051, 5264462.701132784, 2019.0)], | |
203 | ) | |
204 | ||
205 | ||
206 | def test_2d_with_time_itransform_original_crs_obs2(): | |
207 | assert_almost_equal( | |
208 | list( | |
209 | itransform(4896, 7930, [(3496737.2679, 743254.4507, 2019.0)], time_3rd=True) | |
210 | ), | |
211 | [(3496737.4105305015, 743254.1014318303, 2019.0)], | |
212 | ) | |
213 | ||
214 | ||
215 | def test_itransform_time_3rd_invalid(): | |
216 | ||
217 | with pytest.raises(ValueError, match="'time_3rd' is only valid for 3 coordinates."): | |
218 | list( | |
219 | itransform( | |
220 | 7789, | |
221 | 8401, | |
222 | [(3496737.2679, 743254.4507, 5264462.9620, 2019.0)], | |
223 | time_3rd=True, | |
224 | ) | |
225 | ) | |
226 | with pytest.raises(ValueError, match="'time_3rd' is only valid for 3 coordinates."): | |
227 | list(itransform(7789, 8401, [(3496737.2679, 743254.4507)], time_3rd=True)) | |
228 | ||
229 | ||
230 | def test_transform_error(): | |
231 | pj = Proj(init="epsg:4555") | |
232 | pjx, pjy = pj(116.366, 39.867) | |
233 | with pytest.raises(ProjError): | |
234 | transform(pj, Proj(4326), pjx, pjy, radians=True, errcheck=True) | |
235 | ||
236 | ||
237 | def test_itransform_error(): | |
238 | pj = Proj(init="epsg:4555") | |
239 | pjx, pjy = pj(116.366, 39.867) | |
240 | with pytest.raises(ProjError): | |
241 | list(itransform(pj, Proj(4326), [(pjx, pjy)], radians=True, errcheck=True)) | |
242 | ||
243 | ||
244 | def test_transform_radians(): | |
245 | WGS84 = pyproj.Proj("+init=EPSG:4326") | |
246 | ECEF = pyproj.Proj(proj="geocent", ellps="WGS84", datum="WGS84") | |
247 | assert_almost_equal( | |
248 | pyproj.transform( | |
249 | ECEF, WGS84, -2704026.010, -4253051.810, 3895878.820, radians=True | |
250 | ), | |
251 | (-2.137113493845668, 0.6613203738996222, -20.531156923621893), | |
252 | ) | |
253 | ||
254 | assert_almost_equal( | |
255 | pyproj.transform( | |
256 | WGS84, | |
257 | ECEF, | |
258 | -2.137113493845668, | |
259 | 0.6613203738996222, | |
260 | -20.531156923621893, | |
261 | radians=True, | |
262 | ), | |
263 | (-2704026.010, -4253051.810, 3895878.820), | |
264 | ) | |
265 | ||
266 | ||
267 | def test_itransform_radians(): | |
268 | WGS84 = pyproj.Proj("+init=EPSG:4326") | |
269 | ECEF = pyproj.Proj(proj="geocent", ellps="WGS84", datum="WGS84") | |
270 | assert_almost_equal( | |
271 | list( | |
272 | pyproj.itransform( | |
273 | ECEF, WGS84, [(-2704026.010, -4253051.810, 3895878.820)], radians=True | |
274 | ) | |
275 | ), | |
276 | [(-2.137113493845668, 0.6613203738996222, -20.531156923621893)], | |
277 | ) | |
278 | ||
279 | assert_almost_equal( | |
280 | list( | |
281 | pyproj.itransform( | |
282 | WGS84, | |
283 | ECEF, | |
284 | [(-2.137113493845668, 0.6613203738996222, -20.531156923621893)], | |
285 | radians=True, | |
286 | ) | |
287 | ), | |
288 | [(-2704026.010, -4253051.810, 3895878.820)], | |
289 | ) |
0 | # -*- coding: utf-8 -*- | |
1 | """Rewrite part of test.py in pyproj in the form of unittests.""" | |
2 | ||
3 | import math | |
4 | import unittest | |
5 | ||
6 | from pyproj import Geod, Proj, pj_ellps, pj_list, transform | |
7 | from pyproj.crs import CRSError | |
8 | ||
9 | try: | |
10 | import nose2 | |
11 | import nose2.tools | |
12 | ||
13 | HAS_NOSE2 = True | |
14 | except ImportError: | |
15 | HAS_NOSE2 = False | |
16 | ||
17 | ||
18 | class BasicTest(unittest.TestCase): | |
19 | def testProj4Version(self): | |
20 | awips221 = Proj(proj="lcc", R=6371200, lat_1=50, lat_2=50, lon_0=-107) | |
21 | # self.assertEqual(awips221.proj_version, 4.9) | |
22 | ||
23 | def testInitWithBackupString4(self): | |
24 | # this fails unless backup of to_string(4) is used | |
25 | pj = Proj( | |
26 | "+proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@null +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs" | |
27 | ) | |
28 | assert pj.crs.is_valid | |
29 | ||
30 | def testProjAwips221(self): | |
31 | # AWIPS is Advanced Weather Interactive Processing System | |
32 | params = {"proj": "lcc", "R": 6371200, "lat_1": 50, "lat_2": 50, "lon_0": -107} | |
33 | nx = 349 | |
34 | ny = 277 | |
35 | awips221 = Proj( | |
36 | proj=params["proj"], | |
37 | R=params["R"], | |
38 | lat_1=params["lat_1"], | |
39 | lat_2=params["lat_2"], | |
40 | lon_0=params["lon_0"], | |
41 | preserve_units=False, | |
42 | ) | |
43 | awips221_from_dict = Proj(params, preserve_units=False) | |
44 | ||
45 | items = sorted([val for val in awips221.crs.srs.split() if val]) | |
46 | items_dict = sorted([val for val in awips221_from_dict.crs.srs.split() if val]) | |
47 | self.assertEqual(items, items_dict) | |
48 | ||
49 | expected = sorted( | |
50 | [ | |
51 | "+proj=lcc", | |
52 | "+R=6371200", | |
53 | "+lat_1=50", | |
54 | "+lat_2=50", | |
55 | "+lon_0=-107", | |
56 | "+type=crs", | |
57 | ] | |
58 | ) | |
59 | self.assertEqual(items, expected) | |
60 | ||
61 | point = awips221(-145.5, 1.0) | |
62 | x, y = -5632642.22547495, 1636571.4883145525 | |
63 | self.assertAlmostEqual(point[0], x) | |
64 | self.assertAlmostEqual(point[1], y) | |
65 | ||
66 | pairs = [ | |
67 | [(-45, 45), (4351601.20766915, 7606948.029327129)], | |
68 | [(45, 45), (5285389.07739382, 14223336.17467613)], | |
69 | [(45, -45), (20394982.466924712, 21736546.456803113)], | |
70 | [(-45, -45), (16791730.756976362, -3794425.4816524936)], | |
71 | ] | |
72 | for point_geog, expected in pairs: | |
73 | point = awips221(*point_geog) | |
74 | self.assertAlmostEqual(point[0], expected[0]) | |
75 | self.assertAlmostEqual(point[1], expected[1]) | |
76 | point_geog2 = awips221(*point, inverse=True) | |
77 | self.assertAlmostEqual(point_geog[0], point_geog2[0]) | |
78 | self.assertAlmostEqual(point_geog[1], point_geog2[1]) | |
79 | ||
80 | def test_from_dict_with_bool(self): | |
81 | # issue #183 | |
82 | p_d = {'proj': 'omerc', | |
83 | 'lat_2': 80.27942, | |
84 | 'lat_0': 62.87671, | |
85 | 'lat_1': 42.751232, | |
86 | 'ellps': 'WGS84', | |
87 | 'no_rot': True, | |
88 | 'lon_1': 33.793186, | |
89 | 'lon_2': -18.374414} | |
90 | p=Proj(p_d) | |
91 | self.assertTrue('+no_rot' in p.srs.split()) | |
92 | p_d = {'proj': 'omerc', | |
93 | 'lat_2': 80.27942, | |
94 | 'lat_0': 62.87671, | |
95 | 'lat_1': 42.751232, | |
96 | 'ellps': 'WGS84', | |
97 | 'no_rot': False, | |
98 | 'lon_1': 33.793186, | |
99 | 'lon_2': -18.374414} | |
100 | p=Proj(p_d) | |
101 | self.assertFalse('+no_rot' in p.srs.split()) | |
102 | ||
103 | ||
104 | ||
105 | class InverseHammerTest(unittest.TestCase): | |
106 | # This is a unit test of the inverse of the hammer projection, which | |
107 | # was added in the 4.9.3 version of PROJ (then PROJ.4). | |
108 | @classmethod | |
109 | def setUpClass(self): | |
110 | self.p = Proj(proj="hammer") # hammer proj | |
111 | self.x, self.y = self.p(-30, 40) | |
112 | ||
113 | def test_forward(self): | |
114 | self.assertAlmostEqual(self.x, -2711575.083, places=3) | |
115 | self.assertAlmostEqual(self.y, 4395506.619, places=3) | |
116 | ||
117 | def test_inverse(self): | |
118 | lon, lat = self.p(self.x, self.y, inverse=True) | |
119 | self.assertAlmostEqual(lon, -30.0, places=3) | |
120 | self.assertAlmostEqual(lat, 40.0, places=3) | |
121 | ||
122 | ||
123 | class TypeError_Transform_Issue8_Test(unittest.TestCase): | |
124 | # Test for "Segmentation fault on pyproj.transform #8" | |
125 | # https://github.com/jswhit/pyproj/issues/8 | |
126 | ||
127 | def setUp(self): | |
128 | self.p = Proj(init="epsg:4269") | |
129 | ||
130 | def test_tranform_none_1st_parmeter(self): | |
131 | # test should raise Type error if projections are not of Proj classes | |
132 | # version 1.9.4 produced AttributeError, now should raise TypeError | |
133 | with self.assertRaises(CRSError): | |
134 | transform(None, self.p, -74, 39) | |
135 | ||
136 | def test_tranform_none_2nd_parmeter(self): | |
137 | # test should raise Type error if projections are not of Proj classes | |
138 | # version 1.9.4 has a Segmentation Fault, now should raise TypeError | |
139 | with self.assertRaises(CRSError): | |
140 | transform(self.p, None, -74, 39) | |
141 | ||
142 | ||
143 | class Geod_NoDefs_Issue22_Test(unittest.TestCase): | |
144 | # Test for Issue #22, Geod with "+no_defs" in initstring | |
145 | # Before PR #23 merged 2015-10-07, having +no_defs in the initstring would result in a ValueError | |
146 | def test_geod_nodefs(self): | |
147 | Geod("+a=6378137 +b=6378137 +no_defs") | |
148 | ||
149 | ||
150 | class ProjLatLongTypeErrorTest(unittest.TestCase): | |
151 | # .latlong() using in transform raised a TypeError in release 1.9.5.1 | |
152 | # reported in issue #53, resolved in #73. | |
153 | def test_latlong_typeerror(self): | |
154 | p = Proj("+proj=stere +lon_0=-39 +lat_0=90 +lat_ts=71.0 +ellps=WGS84") | |
155 | self.assertTrue(isinstance(p, Proj)) | |
156 | # if not patched this line raises a "TypeError: p2 must be a Proj class" | |
157 | lon, lat = transform(p, p.to_latlong(), 200000, 400000) | |
158 | ||
159 | ||
160 | @unittest.skipUnless(HAS_NOSE2, "nose2 is not installed") | |
161 | class ForwardInverseTest(unittest.TestCase): | |
162 | @nose2.tools.params(*pj_list.keys()) | |
163 | def test_fwd_inv(self, pj): | |
164 | try: | |
165 | p = Proj(proj=pj) | |
166 | x, y = p(-30, 40) | |
167 | # note, for proj 4.9.2 or before the inverse projection may be missing | |
168 | # and pyproj 1.9.5.1 or before does not test for this and will | |
169 | # give a segmentation fault at this point: | |
170 | lon, lat = p(x, y, inverse=True) | |
171 | except RuntimeError: | |
172 | pass | |
173 | ||
174 | ||
175 | # Tests for shared memory between Geod objects | |
176 | class GeodSharedMemoryBugTestIssue64(unittest.TestCase): | |
177 | def setUp(self): | |
178 | self.g = Geod(ellps="clrk66") | |
179 | self.ga = self.g.a | |
180 | self.mercury = Geod(a=2439700) # Mercury 2000 ellipsoid | |
181 | # Mercury is much smaller than earth. | |
182 | ||
183 | def test_not_shared_memory(self): | |
184 | self.assertEqual(self.ga, self.g.a) | |
185 | # mecury must have a different major axis from earth | |
186 | self.assertNotEqual(self.g.a, self.mercury.a) | |
187 | self.assertNotEqual(self.g.b, self.mercury.b) | |
188 | self.assertNotEqual(self.g.sphere, self.mercury.sphere) | |
189 | self.assertNotEqual(self.g.f, self.mercury.f) | |
190 | self.assertNotEqual(self.g.es, self.mercury.es) | |
191 | ||
192 | # initstrings were not shared in issue #64 | |
193 | self.assertNotEqual(self.g.initstring, self.mercury.initstring) | |
194 | ||
195 | def test_distances(self): | |
196 | # note calculated distance was not an issue with #64, but it still a shared memory test | |
197 | boston_lat = 42.0 + (15.0 / 60.0) | |
198 | boston_lon = -71.0 - (7.0 / 60.0) | |
199 | portland_lat = 45.0 + (31.0 / 60.0) | |
200 | portland_lon = -123.0 - (41.0 / 60.0) | |
201 | ||
202 | az12, az21, dist_g = self.g.inv( | |
203 | boston_lon, boston_lat, portland_lon, portland_lat | |
204 | ) | |
205 | ||
206 | az12, az21, dist_mercury = self.mercury.inv( | |
207 | boston_lon, boston_lat, portland_lon, portland_lat | |
208 | ) | |
209 | self.assertLess(dist_mercury, dist_g) | |
210 | ||
211 | ||
212 | class ReprTests(unittest.TestCase): | |
213 | # test __repr__ for Proj object | |
214 | def test_repr(self): | |
215 | p = Proj(proj="latlong", preserve_units=True) | |
216 | expected = "Proj('+proj=longlat +datum=WGS84 +no_defs', preserve_units=True)" | |
217 | self.assertEqual(repr(p), expected) | |
218 | ||
219 | # test __repr__ for Geod object | |
220 | def test_sphere(self): | |
221 | # ellipse is Venus 2000 (IAU2000:29900), which is a sphere | |
222 | g = Geod("+a=6051800 +b=6051800") | |
223 | self.assertEqual(repr(g), "Geod('+a=6051800 +f=0')") | |
224 | ||
225 | # test __repr__ for Geod object | |
226 | def test_ellps_name_round_trip(self): | |
227 | # this could be done in a parameter fashion | |
228 | for ellps_name in pj_ellps: | |
229 | # skip tests, these ellipses NWL9D and WGS66 are the same | |
230 | if ellps_name in ("NWL9D", "WGS66"): | |
231 | continue | |
232 | p = Geod(ellps=ellps_name) | |
233 | expected = "Geod(ellps='{0}')".format(ellps_name) | |
234 | self.assertEqual(repr(p), expected) | |
235 | ||
236 | ||
237 | class TestRadians(unittest.TestCase): | |
238 | """Tests issue #84""" | |
239 | ||
240 | def setUp(self): | |
241 | self.g = Geod(ellps="clrk66") | |
242 | self.boston_d = (-71.0 - (7.0 / 60.0), 42.0 + (15.0 / 60.0)) | |
243 | self.boston_r = (math.radians(self.boston_d[0]), math.radians(self.boston_d[1])) | |
244 | self.portland_d = (-123.0 - (41.0 / 60.0), 45.0 + (31.0 / 60.0)) | |
245 | self.portland_r = ( | |
246 | math.radians(self.portland_d[0]), | |
247 | math.radians(self.portland_d[1]), | |
248 | ) | |
249 | ||
250 | def test_inv_radians(self): | |
251 | ||
252 | # Get bearings and distance from Boston to Portland in degrees | |
253 | az12_d, az21_d, dist_d = self.g.inv( | |
254 | self.boston_d[0], | |
255 | self.boston_d[1], | |
256 | self.portland_d[0], | |
257 | self.portland_d[1], | |
258 | radians=False, | |
259 | ) | |
260 | ||
261 | # Get bearings and distance from Boston to Portland in radians | |
262 | az12_r, az21_r, dist_r = self.g.inv( | |
263 | self.boston_r[0], | |
264 | self.boston_r[1], | |
265 | self.portland_r[0], | |
266 | self.portland_r[1], | |
267 | radians=True, | |
268 | ) | |
269 | ||
270 | # Check they are equal | |
271 | self.assertAlmostEqual(az12_d, math.degrees(az12_r)) | |
272 | self.assertAlmostEqual(az21_d, math.degrees(az21_r)) | |
273 | self.assertAlmostEqual(dist_d, dist_r) | |
274 | ||
275 | def test_fwd_radians(self): | |
276 | # Get bearing and distance to Portland | |
277 | az12_d, az21_d, dist = self.g.inv( | |
278 | self.boston_d[0], | |
279 | self.boston_d[1], | |
280 | self.portland_d[0], | |
281 | self.portland_d[1], | |
282 | radians=False, | |
283 | ) | |
284 | ||
285 | # Calculate Portland's lon/lat from bearing and distance in degrees | |
286 | endlon_d, endlat_d, backaz_d = self.g.fwd( | |
287 | self.boston_d[0], self.boston_d[1], az12_d, dist, radians=False | |
288 | ) | |
289 | ||
290 | # Calculate Portland's lon/lat from bearing and distance in radians | |
291 | endlon_r, endlat_r, backaz_r = self.g.fwd( | |
292 | self.boston_r[0], self.boston_r[1], math.radians(az12_d), dist, radians=True | |
293 | ) | |
294 | ||
295 | # Check they are equal | |
296 | self.assertAlmostEqual(endlon_d, math.degrees(endlon_r)) | |
297 | self.assertAlmostEqual(endlat_d, math.degrees(endlat_r)) | |
298 | self.assertAlmostEqual(backaz_d, math.degrees(backaz_r)) | |
299 | ||
300 | # Check to make sure we're back in Portland | |
301 | self.assertAlmostEqual(endlon_d, self.portland_d[0]) | |
302 | self.assertAlmostEqual(endlat_d, self.portland_d[1]) | |
303 | ||
304 | def test_npts_radians(self): | |
305 | # Calculate 10 points between Boston and Portland in degrees | |
306 | points_d = self.g.npts( | |
307 | lon1=self.boston_d[0], | |
308 | lat1=self.boston_d[1], | |
309 | lon2=self.portland_d[0], | |
310 | lat2=self.portland_d[1], | |
311 | npts=10, | |
312 | radians=False, | |
313 | ) | |
314 | ||
315 | # Calculate 10 points between Boston and Portland in radians | |
316 | points_r = self.g.npts( | |
317 | lon1=self.boston_r[0], | |
318 | lat1=self.boston_r[1], | |
319 | lon2=self.portland_r[0], | |
320 | lat2=self.portland_r[1], | |
321 | npts=10, | |
322 | radians=True, | |
323 | ) | |
324 | ||
325 | # Check they are equal | |
326 | for index, dpoint in enumerate(points_d): | |
327 | self.assertAlmostEqual(dpoint[0], math.degrees(points_r[index][0])) | |
328 | self.assertAlmostEqual(dpoint[1], math.degrees(points_r[index][1])) | |
329 | ||
330 | ||
331 | class Geod_NaN_Issue112_Test(unittest.TestCase): | |
332 | # Test for Issue #112; Geod should silently propagate NaNs in input | |
333 | # to the output. | |
334 | def test_geod_nans(self): | |
335 | g = Geod(ellps="clrk66") | |
336 | (azi1, azi2, s12) = g.inv(43, 10, float("nan"), 20) | |
337 | self.assertTrue(azi1 != azi1) | |
338 | self.assertTrue(azi2 != azi2) | |
339 | self.assertTrue(s12 != s12) | |
340 | (azi1, azi2, s12) = g.inv(43, 10, 53, float("nan")) | |
341 | self.assertTrue(azi1 != azi1) | |
342 | self.assertTrue(azi2 != azi2) | |
343 | self.assertTrue(s12 != s12) | |
344 | # Illegal latitude is treated as NaN | |
345 | (azi1, azi2, s12) = g.inv(43, 10, 53, 91) | |
346 | self.assertTrue(azi1 != azi1) | |
347 | self.assertTrue(azi2 != azi2) | |
348 | self.assertTrue(s12 != s12) | |
349 | (lon2, lat2, azi2) = g.fwd(43, 10, float("nan"), 1e6) | |
350 | self.assertTrue(lon2 != lon2) | |
351 | self.assertTrue(lat2 != lat2) | |
352 | self.assertTrue(azi2 != azi2) | |
353 | (lon2, lat2, azi2) = g.fwd(43, 10, 20, float("nan")) | |
354 | self.assertTrue(lon2 != lon2) | |
355 | self.assertTrue(lat2 != lat2) | |
356 | self.assertTrue(azi2 != azi2) | |
357 | (lon2, lat2, azi2) = g.fwd(43, float("nan"), 20, 1e6) | |
358 | self.assertTrue(lon2 != lon2) | |
359 | self.assertTrue(lat2 != lat2) | |
360 | self.assertTrue(azi2 != azi2) | |
361 | # Illegal latitude is treated as NaN | |
362 | (lon2, lat2, azi2) = g.fwd(43, 91, 20, 1e6) | |
363 | self.assertTrue(lon2 != lon2) | |
364 | self.assertTrue(lat2 != lat2) | |
365 | self.assertTrue(azi2 != azi2) | |
366 | # Only lon2 is NaN | |
367 | (lon2, lat2, azi2) = g.fwd(float("nan"), 10, 20, 1e6) | |
368 | self.assertTrue(lon2 != lon2) | |
369 | self.assertTrue(lat2 == lat2) | |
370 | self.assertTrue(azi2 == azi2) | |
371 | ||
372 | ||
373 | def test_proj_equals(): | |
374 | assert Proj(4326) == Proj("epsg:4326") | |
375 | assert Proj(4326) != Proj("epsg:3857") | |
376 | assert Proj(4326) == Proj(Proj("epsg:4326").crs.to_proj4()) | |
377 | ||
378 | ||
379 | if __name__ == "__main__": | |
380 | if HAS_NOSE2 is True: | |
381 | nose2.discover() | |
382 | else: | |
383 | unittest.main() |
0 | import pytest | |
1 | ||
2 | from pyproj import CRS | |
3 | from pyproj.exceptions import CRSError | |
4 | ||
5 | ||
6 | def test_from_proj4_json(): | |
7 | json_str = '{"proj": "longlat", "ellps": "WGS84", "datum": "WGS84"}' | |
8 | proj = CRS.from_string(json_str) | |
9 | assert proj.to_proj4(4) == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
10 | assert proj.to_proj4(5) == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
11 | # Test with invalid JSON code | |
12 | with pytest.raises(CRSError): | |
13 | assert CRS.from_string("{foo: bar}") | |
14 | ||
15 | ||
16 | def test_from_epsg(): | |
17 | proj = CRS.from_epsg(4326) | |
18 | assert proj.to_epsg() == 4326 | |
19 | ||
20 | # Test with invalid EPSG code | |
21 | with pytest.raises(CRSError): | |
22 | assert CRS.from_epsg(0) | |
23 | ||
24 | ||
25 | def test_from_epsg_string(): | |
26 | proj = CRS.from_string("epsg:4326") | |
27 | assert proj.to_epsg() == 4326 | |
28 | ||
29 | # Test with invalid EPSG code | |
30 | with pytest.raises(ValueError): | |
31 | assert CRS.from_string("epsg:xyz") | |
32 | ||
33 | ||
34 | def test_from_string(): | |
35 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
36 | assert wgs84_crs.to_proj4() == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
37 | # Make sure this doesn't get handled using the from_epsg() even though 'epsg' is in the string | |
38 | epsg_init_crs = CRS.from_string("+init=epsg:26911 +units=m +no_defs=True") | |
39 | assert ( | |
40 | epsg_init_crs.to_proj4() | |
41 | == "+proj=utm +zone=11 +datum=NAD83 +units=m +no_defs +type=crs" | |
42 | ) | |
43 | ||
44 | ||
45 | def test_bare_parameters(): | |
46 | """ Make sure that bare parameters (e.g., no_defs) are handled properly, | |
47 | even if they come in with key=True. This covers interaction with pyproj, | |
48 | which makes presents bare parameters as key=<bool>.""" | |
49 | ||
50 | # Example produced by pyproj | |
51 | proj = CRS.from_string( | |
52 | "+proj=lcc +lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
53 | ) | |
54 | assert "+no_defs" in proj.to_proj4(4) | |
55 | ||
56 | # TODO: THIS DOES NOT WORK | |
57 | proj = CRS.from_string( | |
58 | "+lon_0=-95 +ellps=GRS80 +proj=lcc +y_0=0 +no_defs=False +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
59 | ) | |
60 | # assert "+no_defs" not in proj.to_proj4(4) | |
61 | ||
62 | ||
63 | def test_is_geographic(): | |
64 | assert CRS({"init": "EPSG:4326"}).is_geographic is True | |
65 | assert CRS({"init": "EPSG:3857"}).is_geographic is False | |
66 | ||
67 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
68 | assert wgs84_crs.is_geographic is True | |
69 | ||
70 | nad27_crs = CRS.from_string("+proj=longlat +ellps=clrk66 +datum=NAD27 +no_defs") | |
71 | assert nad27_crs.is_geographic is True | |
72 | ||
73 | lcc_crs = CRS.from_string( | |
74 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
75 | ) | |
76 | assert lcc_crs.is_geographic is False | |
77 | ||
78 | ||
79 | def test_is_projected(): | |
80 | assert CRS({"init": "EPSG:3857"}).is_projected is True | |
81 | ||
82 | lcc_crs = CRS.from_string( | |
83 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
84 | ) | |
85 | assert CRS.from_user_input(lcc_crs).is_projected is True | |
86 | ||
87 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
88 | assert CRS.from_user_input(wgs84_crs).is_projected is False | |
89 | ||
90 | ||
91 | def test_is_same_crs(): | |
92 | crs1 = CRS({"init": "EPSG:4326"}) | |
93 | crs2 = CRS({"init": "EPSG:3857"}) | |
94 | ||
95 | assert crs1 == crs1 | |
96 | assert crs1 != crs2 | |
97 | ||
98 | wgs84_crs = CRS.from_string("+proj=longlat +ellps=WGS84 +datum=WGS84") | |
99 | assert crs1 == wgs84_crs | |
100 | ||
101 | # Make sure that same projection with different parameter are not equal | |
102 | lcc_crs1 = CRS.from_string( | |
103 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
104 | ) | |
105 | lcc_crs2 = CRS.from_string( | |
106 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc +x_0=0 +units=m +lat_2=77 +lat_1=45 +lat_0=0" | |
107 | ) | |
108 | assert lcc_crs1 != lcc_crs2 | |
109 | ||
110 | ||
111 | def test_to_proj4(): | |
112 | assert ( | |
113 | CRS({"init": "EPSG:4326"}).to_proj4(4) | |
114 | == "+proj=longlat +datum=WGS84 +no_defs +type=crs" | |
115 | ) | |
116 | ||
117 | ||
118 | def test_is_valid_false(): | |
119 | with pytest.raises(CRSError): | |
120 | CRS(init="EPSG:432600").is_valid | |
121 | ||
122 | ||
123 | def test_is_valid(): | |
124 | assert CRS(init="EPSG:4326").is_valid | |
125 | ||
126 | ||
127 | def test_empty_json(): | |
128 | with pytest.raises(CRSError): | |
129 | CRS.from_string("{}") | |
130 | with pytest.raises(CRSError): | |
131 | CRS.from_string("[]") | |
132 | with pytest.raises(CRSError): | |
133 | CRS.from_string("") | |
134 | ||
135 | ||
136 | def test_has_wkt_property(): | |
137 | assert ( | |
138 | CRS({"init": "EPSG:4326"}) | |
139 | .to_wkt("WKT1_GDAL") | |
140 | .startswith('GEOGCS["WGS 84",DATUM') | |
141 | ) | |
142 | ||
143 | ||
144 | def test_repr(): | |
145 | assert repr(CRS({"init": "EPSG:4326"})) == ( | |
146 | "<CRS: +init=epsg:4326 +type=crs>\n" | |
147 | "Name: WGS 84\n" | |
148 | "Ellipsoid:\n" | |
149 | "- semi_major_metre: 6378137.00\n" | |
150 | "- semi_minor_metre: 6356752.31\n" | |
151 | "- inverse_flattening: 298.26\n" | |
152 | "Area of Use:\n" | |
153 | "- name: World\n" | |
154 | "- bounds: (-180.0, -90.0, 180.0, 90.0)\n" | |
155 | "Prime Meridian:\n" | |
156 | "- longitude: 0.0000\n" | |
157 | "- unit_name: degree\n" | |
158 | "- unit_conversion_factor: 0.01745329\n" | |
159 | "Axis Info:\n" | |
160 | "- Longitude[lon] (east) EPSG:9122 (degree)\n" | |
161 | "- Latitude[lat] (north) EPSG:9122 (degree)\n" | |
162 | ) | |
163 | ||
164 | ||
165 | def test_repr__long(): | |
166 | assert repr(CRS(CRS({"init": "EPSG:4326"}).to_wkt())) == ( | |
167 | '<CRS: GEOGCRS["WGS 84",DATUM["World Geodetic System 1984 ...>\n' | |
168 | "Name: WGS 84\n" | |
169 | "Ellipsoid:\n" | |
170 | "- semi_major_metre: 6378137.00\n" | |
171 | "- semi_minor_metre: 6356752.31\n" | |
172 | "- inverse_flattening: 298.26\n" | |
173 | "Area of Use:\n" | |
174 | "- name: World\n" | |
175 | "- bounds: (-180.0, -90.0, 180.0, 90.0)\n" | |
176 | "Prime Meridian:\n" | |
177 | "- longitude: 0.0000\n" | |
178 | "- unit_name: degree\n" | |
179 | "- unit_conversion_factor: 0.01745329\n" | |
180 | "Axis Info:\n" | |
181 | "- Longitude[lon] (east) EPSG:9122 (degree)\n" | |
182 | "- Latitude[lat] (north) EPSG:9122 (degree)\n" | |
183 | ) | |
184 | ||
185 | ||
186 | def test_repr__undefined(): | |
187 | assert repr( | |
188 | CRS( | |
189 | "+proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@null" | |
190 | " +lon_0=0.0 +x_0=0.0 +y_0=0.0 +units=m +no_defs" | |
191 | ) | |
192 | ) == ( | |
193 | "<CRS: +proj=merc +a=6378137.0 +b=6378137.0 +nadgrids=@nu ...>\n" | |
194 | "Name: unknown\n" | |
195 | "Ellipsoid:\n" | |
196 | "- semi_major_metre: 6378137.00\n" | |
197 | "- semi_minor_metre: nan\n" | |
198 | "- inverse_flattening: 0.00\n" | |
199 | "Area of Use:\n" | |
200 | "- UNDEFINED\n" | |
201 | "Prime Meridian:\n" | |
202 | "- longitude: 0.0000\n" | |
203 | "- unit_name: degree\n" | |
204 | "- unit_conversion_factor: 0.01745329\n" | |
205 | "Axis Info:\n" | |
206 | "- UNDEFINED" | |
207 | ) | |
208 | ||
209 | ||
210 | def test_dunder_str(): | |
211 | assert str(CRS({"init": "EPSG:4326"})) == CRS({"init": "EPSG:4326"}).srs | |
212 | ||
213 | ||
214 | def test_epsg(): | |
215 | assert CRS({"init": "EPSG:4326"}).to_epsg(20) == 4326 | |
216 | assert CRS({"init": "EPSG:4326"}).to_epsg() is None | |
217 | assert CRS.from_user_input(4326).to_epsg() == 4326 | |
218 | assert CRS.from_epsg(4326).to_epsg() == 4326 | |
219 | assert CRS.from_user_input("epsg:4326").to_epsg() == 4326 | |
220 | ||
221 | ||
222 | def test_datum(): | |
223 | assert CRS.from_epsg(4326).datum == CRS( | |
224 | 'DATUM["World Geodetic System 1984",' | |
225 | 'ELLIPSOID["WGS 84",6378137,298.257223563,' | |
226 | 'LENGTHUNIT["metre",1]],ID["EPSG",6326]]' | |
227 | ) | |
228 | ||
229 | ||
230 | def test_epsg__not_found(): | |
231 | assert CRS("+proj=longlat +datum=WGS84 +no_defs").to_epsg(0) is None | |
232 | assert CRS.from_string("+proj=longlat +datum=WGS84 +no_defs").to_epsg() is None | |
233 | ||
234 | ||
235 | def test_epsg__no_code_available(): | |
236 | lcc_crs = CRS.from_string( | |
237 | "+lon_0=-95 +ellps=GRS80 +y_0=0 +no_defs=True +proj=lcc " | |
238 | "+x_0=0 +units=m +lat_2=77 +lat_1=49 +lat_0=0" | |
239 | ) | |
240 | assert lcc_crs.to_epsg() is None | |
241 | ||
242 | ||
243 | def test_crs_OSR_equivalence(): | |
244 | crs1 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
245 | crs2 = CRS.from_string("+proj=latlong +datum=WGS84 +no_defs") | |
246 | crs3 = CRS({"init": "EPSG:4326"}) | |
247 | assert crs1 == crs2 | |
248 | # these are not equivalent in proj.4 now as one uses degrees and the othe radians | |
249 | assert crs1 == crs3 | |
250 | ||
251 | ||
252 | def test_crs_OSR_no_equivalence(): | |
253 | crs1 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
254 | crs2 = CRS.from_string("+proj=longlat +datum=NAD27 +no_defs") | |
255 | assert crs1 != crs2 | |
256 | ||
257 | ||
258 | def test_from_wkt(): | |
259 | wgs84 = CRS.from_string("+proj=longlat +datum=WGS84 +no_defs") | |
260 | from_wkt = CRS(wgs84.to_wkt()) | |
261 | assert wgs84.to_wkt() == from_wkt.to_wkt() | |
262 | ||
263 | ||
264 | def test_from_wkt_invalid(): | |
265 | with pytest.raises(CRSError): | |
266 | CRS("trash-54322") | |
267 | with pytest.raises(CRSError): | |
268 | CRS("") | |
269 | ||
270 | ||
271 | def test_from_user_input_epsg(): | |
272 | assert "+proj=longlat" in CRS.from_user_input("EPSG:4326").to_proj4(4) | |
273 | ||
274 | ||
275 | def test_from_esri_wkt(): | |
276 | projection_string = ( | |
277 | 'PROJCS["USA_Contiguous_Albers_Equal_Area_Conic_USGS_version",' | |
278 | 'GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",' | |
279 | 'SPHEROID["GRS_1980",6378137.0,298.257222101]],' | |
280 | 'PRIMEM["Greenwich",0.0],' | |
281 | 'UNIT["Degree",0.0174532925199433]],' | |
282 | 'PROJECTION["Albers"],' | |
283 | 'PARAMETER["false_easting",0.0],' | |
284 | 'PARAMETER["false_northing",0.0],' | |
285 | 'PARAMETER["central_meridian",-96.0],' | |
286 | 'PARAMETER["standard_parallel_1",29.5],' | |
287 | 'PARAMETER["standard_parallel_2",45.5],' | |
288 | 'PARAMETER["latitude_of_origin",23.0],' | |
289 | 'UNIT["Meter",1.0],' | |
290 | 'VERTCS["NAVD_1988",' | |
291 | 'VDATUM["North_American_Vertical_Datum_1988"],' | |
292 | 'PARAMETER["Vertical_Shift",0.0],' | |
293 | 'PARAMETER["Direction",1.0],UNIT["Centimeter",0.01]]]' | |
294 | ) | |
295 | proj_crs_str = CRS.from_string(projection_string) | |
296 | proj_crs_wkt = CRS(projection_string) | |
297 | assert proj_crs_str.to_proj4() == proj_crs_wkt.to_proj4() | |
298 | assert proj_crs_str.to_proj4(4) == ( | |
299 | "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 " | |
300 | "+lat_2=45.5 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs" | |
301 | ) | |
302 | ||
303 | ||
304 | def test_compound_crs(): | |
305 | wkt = """COMPD_CS["unknown",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],VERT_CS["unknown",VERT_DATUM["unknown",2005],UNIT["metre",1.0,AUTHORITY["EPSG","9001"]],AXIS["Up",UP]]]""" | |
306 | assert CRS(wkt).to_wkt("WKT1_GDAL").startswith('COMPD_CS["unknown",GEOGCS["WGS 84"') | |
307 | ||
308 | ||
309 | def test_ellipsoid(): | |
310 | crs1 = CRS.from_epsg(4326) | |
311 | assert "{:.3f}".format(crs1.ellipsoid.inverse_flattening) == "298.257" | |
312 | assert "{:.3f}".format(crs1.ellipsoid.semi_major_metre) == "6378137.000" | |
313 | assert "{:.3f}".format(crs1.ellipsoid.semi_minor_metre) == "6356752.314" | |
314 | ||
315 | ||
316 | def test_area_of_use(): | |
317 | crs1 = CRS.from_epsg(4326) | |
318 | assert crs1.area_of_use.bounds == (-180.0, -90.0, 180.0, 90.0) | |
319 | assert crs1.area_of_use.name == "World" | |
320 | ||
321 | ||
322 | def test_from_user_input_custom_crs_class(): | |
323 | class CustomCRS(object): | |
324 | def to_wkt(self): | |
325 | return CRS.from_epsg(4326).to_wkt() | |
326 | ||
327 | assert CRS.from_user_input(CustomCRS()) == CRS.from_epsg(4326) |
0 | import os | |
1 | import shutil | |
2 | import tempfile | |
3 | import unittest | |
4 | from contextlib import contextmanager | |
5 | ||
6 | import pytest | |
7 | from mock import patch | |
8 | ||
9 | from pyproj.datadir import DataDirError, append_data_dir, get_data_dir, set_data_dir | |
10 | ||
11 | ||
12 | def create_projdb(tmpdir): | |
13 | with open(os.path.join(tmpdir, "proj.db"), "w") as pjdb: | |
14 | pjdb.write("DUMMY proj.db") | |
15 | ||
16 | ||
17 | @contextmanager | |
18 | def proj_env(): | |
19 | """ | |
20 | Ensure environment variable the same at the end of the test. | |
21 | """ | |
22 | proj_lib = os.environ.get("PROJ_LIB") | |
23 | try: | |
24 | yield | |
25 | finally: | |
26 | if proj_lib: | |
27 | os.environ["PROJ_LIB"] = proj_lib | |
28 | ||
29 | ||
30 | @contextmanager | |
31 | def temporary_directory(): | |
32 | """ | |
33 | Get a temporary directory | |
34 | """ | |
35 | temp_dir = tempfile.mkdtemp() | |
36 | try: | |
37 | yield temp_dir | |
38 | finally: | |
39 | shutil.rmtree(temp_dir) | |
40 | ||
41 | ||
42 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
43 | def test_get_data_dir__missing(): | |
44 | with proj_env(), pytest.raises(DataDirError), patch( | |
45 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
46 | ), patch("pyproj.datadir.find_executable", return_value=None): | |
47 | set_data_dir(None) | |
48 | os.environ.pop("PROJ_LIB", None) | |
49 | get_data_dir() | |
50 | ||
51 | ||
52 | def test_get_data_dir__from_user(): | |
53 | with proj_env(), temporary_directory() as tmpdir, temporary_directory() as tmpdir_env: | |
54 | create_projdb(tmpdir) | |
55 | os.environ["PROJ_LIB"] = tmpdir_env | |
56 | create_projdb(tmpdir_env) | |
57 | set_data_dir(tmpdir) | |
58 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
59 | os.makedirs(internal_proj_dir) | |
60 | create_projdb(internal_proj_dir) | |
61 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
62 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
63 | assert get_data_dir() == tmpdir | |
64 | ||
65 | ||
66 | def test_get_data_dir__internal(): | |
67 | with proj_env(), temporary_directory() as tmpdir: | |
68 | set_data_dir(None) | |
69 | os.environ["PROJ_LIB"] = tmpdir | |
70 | create_projdb(tmpdir) | |
71 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
72 | os.makedirs(internal_proj_dir) | |
73 | create_projdb(internal_proj_dir) | |
74 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
75 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
76 | assert get_data_dir() == internal_proj_dir | |
77 | ||
78 | ||
79 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
80 | def test_get_data_dir__from_env_var(): | |
81 | with proj_env(), temporary_directory() as tmpdir, patch( | |
82 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
83 | ): | |
84 | set_data_dir(None) | |
85 | os.environ["PROJ_LIB"] = tmpdir | |
86 | create_projdb(tmpdir) | |
87 | assert get_data_dir() == tmpdir | |
88 | ||
89 | ||
90 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
91 | def test_get_data_dir__from_env_var__multiple(): | |
92 | with proj_env(), temporary_directory() as tmpdir, patch( | |
93 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
94 | ): | |
95 | set_data_dir(None) | |
96 | os.environ["PROJ_LIB"] = ";".join([tmpdir, tmpdir, tmpdir]) | |
97 | create_projdb(tmpdir) | |
98 | assert get_data_dir() == ";".join([tmpdir, tmpdir, tmpdir]) | |
99 | ||
100 | ||
101 | @unittest.skipIf(os.name == "nt", reason="Cannot modify Windows environment variables.") | |
102 | def test_get_data_dir__from_path(): | |
103 | with proj_env(), temporary_directory() as tmpdir, patch( | |
104 | "pyproj.datadir.os.path.abspath", return_value="INVALID" | |
105 | ), patch("pyproj.datadir.find_executable") as find_exe: | |
106 | set_data_dir(None) | |
107 | os.environ.pop("PROJ_LIB", None) | |
108 | find_exe.return_value = os.path.join(tmpdir, "bin", "proj") | |
109 | proj_dir = os.path.join(tmpdir, "share", "proj") | |
110 | os.makedirs(proj_dir) | |
111 | create_projdb(proj_dir) | |
112 | assert get_data_dir() == proj_dir | |
113 | ||
114 | ||
115 | def test_append_data_dir__internal(): | |
116 | with proj_env(), temporary_directory() as tmpdir: | |
117 | set_data_dir(None) | |
118 | os.environ["PROJ_LIB"] = tmpdir | |
119 | create_projdb(tmpdir) | |
120 | internal_proj_dir = os.path.join(tmpdir, "proj_dir", "share", "proj") | |
121 | os.makedirs(internal_proj_dir) | |
122 | create_projdb(internal_proj_dir) | |
123 | extra_datadir = str(os.path.join(tmpdir, "extra_datumgrids")) | |
124 | with patch("pyproj.datadir.os.path.abspath") as abspath_mock: | |
125 | abspath_mock.return_value = os.path.join(tmpdir, "randomfilename.py") | |
126 | append_data_dir(extra_datadir) | |
127 | assert get_data_dir() == ";".join([internal_proj_dir, extra_datadir]) |
0 | """ | |
1 | This is a wrapper for the doctests in lib/pyproj/__init__.py so that | |
2 | nose2 can conveniently run all the tests in a single command line. | |
3 | """ | |
4 | ||
5 | import pyproj | |
6 | ||
7 | ||
8 | def test_doctests(): | |
9 | failure_count = pyproj.test() | |
10 | # if the below line fails, doctests have failed | |
11 | assert ( | |
12 | failure_count == 0 | |
13 | ), "{0} of the doctests in " "lib/pyproj/__init__.py failed".format(failure_count) |
0 | import pytest | |
1 | ||
2 | from pyproj import CRS, Proj | |
3 | from pyproj.exceptions import CRSError, ProjError | |
4 | ||
5 | ||
6 | def test_proj_exception(): | |
7 | with pytest.raises(ProjError, match="Internal Proj Error"): | |
8 | Proj("+proj=bobbyjoe") | |
9 | ||
10 | ||
11 | def test_crs_exception(): | |
12 | with pytest.raises(CRSError, match="Internal Proj Error"): | |
13 | CRS("+proj=bobbyjoe") |
0 | import numpy as np | |
1 | from numpy.testing import assert_almost_equal, assert_equal | |
2 | ||
3 | import pyproj | |
4 | from pyproj import Transformer | |
5 | from pyproj.exceptions import ProjError | |
6 | ||
7 | ||
8 | def test_tranform_wgs84_to_custom(): | |
9 | custom_proj = pyproj.Proj( | |
10 | "+proj=geos +lon_0=0.000000 +lat_0=0 +h=35807.414063" | |
11 | " +a=6378.169000 +b=6356.583984" | |
12 | ) | |
13 | wgs84 = pyproj.Proj("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs") | |
14 | lat, lon = 51.04715, 3.23406 | |
15 | xx, yy = pyproj.transform(wgs84, custom_proj, lon, lat) | |
16 | assert "{:.3f} {:.3f}".format(xx, yy) == "212.623 4604.975" | |
17 | ||
18 | ||
19 | def test_transform_wgs84_to_alaska(): | |
20 | lat_lon_proj = pyproj.Proj(init="epsg:4326", preserve_units=False) | |
21 | alaska_aea_proj = pyproj.Proj(init="epsg:2964", preserve_units=False) | |
22 | test = (-179.72638, 49.752533) | |
23 | xx, yy = pyproj.transform(lat_lon_proj, alaska_aea_proj, *test) | |
24 | assert "{:.3f} {:.3f}".format(xx, yy) == "-1824924.495 330822.800" | |
25 | ||
26 | ||
27 | def test_illegal_transformation(): | |
28 | # issue 202 | |
29 | p1 = pyproj.Proj(init="epsg:4326") | |
30 | p2 = pyproj.Proj(init="epsg:3857") | |
31 | xx, yy = pyproj.transform( | |
32 | p1, p2, (-180, -180, 180, 180, -180), (-90, 90, 90, -90, -90) | |
33 | ) | |
34 | assert np.all(np.isinf(xx)) | |
35 | assert np.all(np.isinf(yy)) | |
36 | try: | |
37 | xx, yy = pyproj.transform( | |
38 | p1, p2, (-180, -180, 180, 180, -180), (-90, 90, 90, -90, -90), errcheck=True | |
39 | ) | |
40 | assert_equal(None, "Should throw an exception when errcheck=True") | |
41 | except ProjError: | |
42 | pass | |
43 | ||
44 | ||
45 | def test_lambert_conformal_transform(): | |
46 | # issue 207 | |
47 | Midelt = pyproj.Proj(init="epsg:26191") | |
48 | WGS84 = pyproj.Proj(init="epsg:4326") | |
49 | ||
50 | E = 567623.931 | |
51 | N = 256422.787 | |
52 | h = 1341.467 | |
53 | ||
54 | Long1, Lat1, H1 = pyproj.transform(Midelt, WGS84, E, N, h, radians=False) | |
55 | assert_almost_equal((Long1, Lat1, H1), (-4.6753456, 32.902199, 1341.467), decimal=5) | |
56 | ||
57 | ||
58 | def test_equivalent_crs(): | |
59 | transformer = Transformer.from_crs("epsg:4326", 4326, skip_equivalent=True) | |
60 | assert transformer._transformer.projections_equivalent | |
61 | assert transformer._transformer.projections_exact_same | |
62 | assert transformer._transformer.skip_equivalent | |
63 | ||
64 | ||
65 | def test_equivalent_crs__disabled(): | |
66 | transformer = Transformer.from_crs("epsg:4326", 4326) | |
67 | assert not transformer._transformer.skip_equivalent | |
68 | assert transformer._transformer.projections_equivalent | |
69 | assert transformer._transformer.projections_exact_same | |
70 | ||
71 | ||
72 | def test_equivalent_crs__different(): | |
73 | transformer = Transformer.from_crs("epsg:4326", 3857, skip_equivalent=True) | |
74 | assert transformer._transformer.skip_equivalent | |
75 | assert not transformer._transformer.projections_equivalent | |
76 | assert not transformer._transformer.projections_exact_same | |
77 | ||
78 | ||
79 | def test_equivalent_proj(): | |
80 | transformer = Transformer.from_proj( | |
81 | "+init=epsg:4326", pyproj.Proj(4326).crs.to_proj4(), skip_equivalent=True | |
82 | ) | |
83 | assert transformer._transformer.skip_equivalent | |
84 | assert transformer._transformer.projections_equivalent | |
85 | assert transformer._transformer.projections_exact_same | |
86 | ||
87 | ||
88 | def test_equivalent_proj__disabled(): | |
89 | transformer = Transformer.from_proj(3857, pyproj.Proj(3857).crs.to_proj4()) | |
90 | assert not transformer._transformer.skip_equivalent | |
91 | assert transformer._transformer.projections_equivalent | |
92 | assert transformer._transformer.projections_exact_same | |
93 | ||
94 | ||
95 | def test_equivalent_proj__different(): | |
96 | transformer = Transformer.from_proj(3857, 4326, skip_equivalent=True) | |
97 | assert transformer._transformer.skip_equivalent | |
98 | assert not transformer._transformer.projections_equivalent | |
99 | assert not transformer._transformer.projections_exact_same | |
100 | ||
101 | ||
102 | def test_equivalent_pipeline(): | |
103 | transformer = Transformer.from_pipeline( | |
104 | "+proj=pipeline +step +proj=longlat +ellps=WGS84 +step " | |
105 | "+proj=unitconvert +xy_in=rad +xy_out=deg" | |
106 | ) | |
107 | assert not transformer._transformer.skip_equivalent | |
108 | assert not transformer._transformer.projections_equivalent | |
109 | assert not transformer._transformer.projections_exact_same |