-
Notifications
You must be signed in to change notification settings - Fork 16
/
geodict.py
executable file
·1056 lines (925 loc) · 33.6 KB
/
geodict.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#!/usr/bin/env python
# third party imports
from curses.ascii import alt
import numpy as np
import pyproj
import rasterio
from rasterio.transform import Affine
from .dataset import DataSetException
GLOBE_SPAN = 359.0
def is_valid_projection(projstr):
try:
pyproj.Proj(projstr)
except BaseException:
return False
return True
class GeoDict(object):
EPS = 1e-12
DEFAULT_PROJ4 = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
DIST_THRESH = 0.01 / (111.191 * 1000) # 1 centimeter in decimal degrees
REQ_KEYS = ["xmin", "xmax", "ymin", "ymax", "dx", "dy", "ny", "nx"]
def __init__(self, geodict, adjust="bounds"):
"""
An object which represents the spatial information for a grid and is guaranteed
to be self-consistent.
:param geodict:
A dictionary containing at least some of the following fields:
- xmin Longitude minimum (decimal degrees) (Center of upper left cell)
- xmax Longitude maximum (decimal degrees) (Center of upper right cell)
- ymin Longitude minimum (decimal degrees) (Center of lower left cell)
- ymax Longitude maximum (decimal degrees) (Center of lower right cell)
- dx Cell width (decimal degrees)
- dy Cell height (decimal degrees)
- ny Number of rows of input data (must match input data dimensions)
- nx Number of columns of input data (must match input data dimensions).
- projection proj4 string defining projection. Optional. If not specified,
assumed to be
"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" (WGS-84 geographic).
:param adjust:
String (one of 'bounds','res')
'bounds': dx/dy, nx/ny, xmin/ymax are assumed to be correct, xmax/ymin
will be recalculated.
'res': nx/ny, xmin/ymax, xmax/ymin and assumed to be correct, dx/dy will
be recalculated.
"""
for key in self.REQ_KEYS:
if key not in geodict.keys():
raise DataSetException(
'Missing required key "%s" from input geodict dictionary' % key
)
self._xmin = float(geodict["xmin"])
self._xmax = float(geodict["xmax"])
self._ymin = float(geodict["ymin"])
self._ymax = float(geodict["ymax"])
self._dx = float(geodict["dx"])
self._dy = float(geodict["dy"])
self._ny = int(geodict["ny"])
self._nx = int(geodict["nx"])
self._nodata = None
if "nodata" in geodict:
self._nodata = geodict["nodata"]
if "projection" in geodict:
projstr = geodict["projection"]
if is_valid_projection(projstr):
self._projection = projstr
else:
raise DataSetException(f"Projection string '{projstr}' is invalid.")
else:
self._projection = self.DEFAULT_PROJ4
self.validate(adjust)
@classmethod
def createDictFromBox(cls, xmin, xmax, ymin, ymax, dx, dy, inside=False):
"""Create GeoDict from two corner points and an x/y resolution.
:param xmin: X coordinate of center of upper left pixel.
:param xmax: X coordinate of center of lower right pixel.
:param ymin: Y coordinate of center of lower right pixel.
:param ymax: Y coordinate of center of upper left pixel.
:param dx: Width of pixels.
:param dy: Height of pixels.
:param inside:
Boolean, indicating whether to ensure that resulting GeoDict
falls inside or outside the input bounds.
:returns:
GeoDict object.
"""
if xmin > xmax:
txmax = xmax + 360
else:
txmax = xmax
if inside:
nx = np.floor(((txmax - xmin + GeoDict.EPS) / dx) + 1)
ny = np.floor(((ymax - ymin + GeoDict.EPS) / dy) + 1)
else:
nx = np.ceil(((txmax - xmin - GeoDict.EPS) / dx) + 1)
ny = np.ceil(((ymax - ymin - GeoDict.EPS) / dy) + 1)
xmax2 = xmin + (nx - 1) * dx
ymin2 = ymax - (ny - 1) * dy
return cls(
{
"xmin": xmin,
"xmax": xmax2,
"ymin": ymin2,
"ymax": ymax,
"dx": dx,
"dy": dy,
"nx": nx,
"ny": ny,
}
)
@classmethod
def createDictFromCenter(cls, cx, cy, dx, dy, xspan, yspan):
"""Create GeoDict from a center point, dx/dy and a width and height.
:param cx: X coordinate of center point.
:param cy: Y coordinate of center point.
:param dx: Width of pixels.
:param dy: Height of pixels.
:param xspan: Width of desired box.
:param yspan: Height of desired box.
:returns:
GeoDict object.
"""
xmin = cx - xspan / 2.0
xmax = cx + xspan / 2.0
ymin = cy - yspan / 2.0
ymax = cy + yspan / 2.0
return cls.createDictFromBox(xmin, xmax, ymin, ymax, dx, dy)
def setProjection(self, projection):
"""Set a new projection for the GeoDict.
:param projection:
Valid proj4 string.
:raises DataSetException:
When input is not valid proj4.
"""
try:
pyproj.Proj(projection)
except BaseException:
raise DataSetException(
"%s is not a valid proj4 string." % geodict["projection"]
)
self._projection = projection
def getAligned(self, geodict, inside=False):
"""Return a geodict that is close to the input geodict bounds but aligned with
this GeoDict.
:param geodict:
Input GeoDict object, whose bounds will be used, but dx/dy and nx/ny ignored.
:param inside:
Boolean indicating whether the aligned geodict should be inside or outside
input geodict.
:returns:
GeoDict object which is guaranteed to be grid-aligned with this GeoDict.
"""
dx = self.dx
dy = self.dy
# how many columns are there between the host and input geodict left edges?
if inside:
falsenx = np.ceil((geodict.xmin - self.xmin) / dx)
else:
falsenx = np.floor((geodict.xmin - self.xmin) / dx)
# use that number of rows to calculate what aligned xmin should be
newxmin = self.xmin + falsenx * dx
# how many columns are there between the host and input geodict right edges?
if inside:
falsenx = np.floor((geodict.xmax - self.xmax) / dx)
else:
falsenx = np.ceil((geodict.xmax - self.xmax) / dx)
# use that number of rows to calculate what aligned xmax should be
newxmax = self.xmax + falsenx * dx
# if we wound up going past 180, correct x values to be within -180->180.
if newxmin > 180:
newxmin -= 360
if newxmax > 180:
newxmax -= 360
# how many columns are there between the host and input geodict bottom edges?
if inside:
falseny = np.ceil((geodict.ymin - self.ymin) / dy)
else:
falseny = np.floor((geodict.ymin - self.ymin) / dy)
# use that number of rows to calculate what aligned ymin should be
newymin = self.ymin + falseny * dy
# how many columns are there between the host and input geodict top edges?
if inside:
falseny = np.floor((geodict.ymax - self.ymax) / dy)
else:
falseny = np.ceil((geodict.ymax - self.ymax) / dy)
# use that number of rows to calculate what aligned ymax should be
newymax = self.ymax + falseny * dy
nx = int(np.round((newxmax - newxmin) / dx + 1))
ny = int(np.round((newymax - newymin) / dy + 1))
gd = GeoDict(
{
"xmin": newxmin,
"xmax": newxmax,
"ymin": newymin,
"ymax": newymax,
"dx": dx,
"dy": dy,
"nx": nx,
"ny": ny,
}
)
return gd
def getIntersection(self, geodict):
"""Return a geodict defining intersected area, retaining resolution of the input
geodict.
:param geodict:
Input GeoDict object, which should intersect with this GeoDict.
:returns:
GeoDict which represents the intersected area, and is aligned with the input
geodict.
:raises DataSetException:
When input geodict does not intersect at all with this GeoDict.
"""
# return modified input geodict where bounds have been adjusted to be inside
# self. output should align with input.
if not self.intersects(geodict):
raise DataSetException("Input geodict has no overlap.")
fxmin, fxmax, fymin, fymax = (self.xmin, self.xmax, self.ymin, self.ymax)
xmin, xmax, ymin, ymax = (
geodict.xmin,
geodict.xmax,
geodict.ymin,
geodict.ymax,
)
# there are some weird special cases here with longitudes (180 meridian)
# handling these in separate testable function
txmin, txmax = get_longitude_intersection(xmin, xmax, fxmin, fxmax)
# calculate latitude intersection
tymin = max(fymin, ymin)
tymax = min(fymax, ymax)
# we're using the input geodict resolution for output
dx, dy = (geodict.dx, geodict.dy)
# now align those bounds with the input geodict
trow, tcol = geodict.getRowCol(tymax, txmin, returnFloat=True)
fleftcol = int(np.round(tcol))
ftoprow = int(np.round(trow))
newymax, newxmin = geodict.getLatLon(ftoprow, fleftcol)
trow, tcol = geodict.getRowCol(tymin, txmax, returnFloat=True)
frightcol = int(np.round(tcol))
fbottomrow = int(np.round(trow))
newymin, newxmax = geodict.getLatLon(fbottomrow, frightcol)
nx = int(np.round((newxmax - newxmin) / dx + 1))
ny = int(np.round((newymax - newymin) / dy + 1))
if newxmax > 180:
newxmax = newxmax - 360
outdict = GeoDict(
{
"xmin": newxmin,
"xmax": newxmax,
"ymin": newymin,
"ymax": newymax,
"dx": dx,
"dy": dy,
"nx": nx,
"ny": ny,
}
)
return outdict
def getBoundsWithin(self, geodict):
"""Create a GeoDict by finding the maximum bounding
box aligned with enclosing GeoDict that is guaranteed
to be inside input GeoDict.
:param geodict:
GeoDict object that output GeoDict will be contained by.
:raises DataSetException:
When input geodict is not fully contained by this GeoDict, or
if the output GeoDict cannot be aligned with this GeoDict (this
shouldn't happen).
"""
if not self.contains(geodict):
raise DataSetException(
"Input geodict not fully contained by this GeoDict object."
)
# if the input geodict is identical to the host grid, then just return
# that geodict.
if self == geodict:
return geodict.copy()
# fxmin etc. are the host bounds
fxmin, fymax = (self.xmin, self.ymax)
# address case where host left edge is positive
if (
self.xmin > geodict.xmin
and self.xmax > geodict.xmin
and self.xmax < self.xmin
):
fxmin -= 360
# xmin etc. are the geodict bounds
xmin, xmax, ymin, ymax = (
geodict.xmin,
geodict.xmax,
geodict.ymin,
geodict.ymax,
)
fdx, fdy = (self.dx, self.dy)
fnx = self.nx
# how many columns from the host xmin is the sample xmin?
xmincol = np.ceil((xmin - fxmin) / fdx)
# how many columns from the host xmin is the sample xmax?
if xmax < fxmin and xmax < 0:
xmaxcol = np.floor(((xmax + 360) - fxmin) / fdx)
else:
xmaxcol = np.floor((xmax - fxmin) / fdx)
newxmin = fxmin + xmincol * fdx
newxmax = fxmin + xmaxcol * fdx
yminrow = np.floor((fymax - ymin) / fdy)
ymaxrow = np.ceil((fymax - ymax) / fdy)
newymin = fymax - yminrow * fdy
newymax = fymax - ymaxrow * fdy
# if our original geodict crossed the 180 meridian, then
# let's set our xmin/xmax back to the xmin > xmax scenario.
if xmax < fxmin:
newxmax -= 360
# in certain edge cases, the "inside" bounds can be equal (or nearly)
# to input bounds. Check for this here, and nudge in one pixel until
# fixed
while newxmin <= xmin:
xmincol = xmincol + 1
newxmin = fxmin + xmincol * fdx
while newxmax >= xmax:
xmaxcol = xmaxcol - 1
newxmax = fxmin + xmaxcol * fdx
while newymin <= ymin:
yminrow = yminrow - 1
newymin = fymax - yminrow * fdy
while newymax >= ymax:
ymaxrow = ymaxrow + 1
newymax = fymax - ymaxrow * fdy
ny = int((yminrow - ymaxrow) + 1)
# if the input geodict crosses the 180 meridian, deal with this...
if xmin > xmax:
nx = ((fnx - xmincol)) + (xmaxcol + 1)
else:
nx = int((xmaxcol - xmincol) + 1)
# because validate() calculates ymin and xmax differently,
# let's re-calculate those that way and see if we get the same results
tmp_new_ymin = newymax - fdy * (ny - 1)
tmp_new_xmax = newxmin + fdx * (nx - 1)
tmp_nx = nx
tmp_ny = ny
tmp_xmax = xmax
if tmp_new_xmax > 180:
tmp_xmax += 360
# if we don't get the same results, adjust nx/ny until we do.
while tmp_new_ymin < ymin:
tmp_ny = tmp_ny - 1
tmp_new_ymin = newymax - fdy * (tmp_ny - 1)
while tmp_new_xmax > tmp_xmax:
tmp_nx = tmp_nx - 1
tmp_new_xmax = newxmin + fdx * (tmp_nx - 1)
if tmp_new_ymin != newymin:
newymin = tmp_new_ymin
ny = tmp_ny
if tmp_new_xmax != newxmax:
newxmax = tmp_new_xmax
nx = tmp_nx
outgeodict = GeoDict(
{
"xmin": newxmin,
"xmax": newxmax,
"ymin": newymin,
"ymax": newymax,
"dx": fdx,
"dy": fdy,
"ny": ny,
"nx": nx,
}
)
isaligned = self.isAligned(outgeodict, altxmin=fxmin)
if not isaligned:
raise DataSetException(
"getBoundsWithin() cannot create an aligned geodict."
)
isinside = geodict.contains(outgeodict)
if not isinside:
raise DataSetException(
"getBoundsWithin() cannot create an geodict inside input."
)
return outgeodict
def isAligned(self, geodict, altxmin=None):
"""Determines if input geodict is grid-aligned with this GeoDict.
:param geodict:
Input geodict whose cells must all be grid-aligned with this GeoDict.
:returns:
True when geodict is grid-aligned, and False if not.
"""
dx1 = self._dx
dx2 = geodict.dx
dy1 = self._dy
dy2 = geodict.dy
dx_close = np.isclose(dx1, dx2, atol=self.EPS)
dy_close = np.isclose(dy1, dy2, atol=self.EPS)
if not dx_close or not dy_close:
return False
xmin1 = self._xmin
if altxmin is not None:
xmin1 = altxmin
xmin2 = geodict.xmin
ymin1 = self._ymin
ymin2 = geodict.ymin
t1 = (xmin2 - xmin1) / dx1
t2 = np.round(t1)
x_close = np.isclose(t1, t2)
t1 = (ymin2 - ymin1) / dy1
t2 = np.round(t1)
y_close = np.isclose(t1, t2)
# x_rem = ((xmin2-xmin1) / dx1) < self.EPS
# y_rem = ((ymin2-ymin1) % dy1) < self.EPS
if x_close and y_close:
return True
return False
def doesNotContain(self, geodict):
"""Determine if input geodict is completely outside this GeoDict.
:param geodict:
Input GeoDict object.
:returns:
True if input geodict is completely outside this GeoDict,
False if not.
"""
a, b = (self.xmin, self.ymax)
e, f = (geodict.xmin, geodict.ymax)
c, _ = (self.xmax, self.ymin)
g, h = (geodict.xmax, geodict.ymin)
outside_x = e > c or a > g
outside_y = f > b or b > h
# outside_x = geodict.xmin < self._xmin and geodict.xmax > self._xmax
# outside_y = geodict.ymin < self._ymin and geodict.ymax > self._ymax
if outside_x and outside_y:
return True
return False
def intersects(self, geodict):
"""Determine if input geodict intersects this GeoDict.
:param geodict:
Input GeoDict object.
:returns:
True if input geodict intersects with this GeoDict,
False if not.
"""
try:
xmin = self.xmin
xmax = self.xmax
fxmin = geodict.xmin
fxmax = geodict.xmax
txmin, txmax = get_longitude_intersection(xmin, xmax, fxmin, fxmax)
except Exception:
return False
ymin = self.ymin
ymax = self.ymax
fymin = geodict.ymin
fymax = geodict.ymax
self_inside_input = ymin <= fymin and ymax >= fymin
input_inside_self = fymin <= ymin and fymax >= ymin
if self_inside_input or input_inside_self:
return True
# if self.xmin > self.xmax:
# c, d = (self.xmax + 360, self.ymin)
# else:
# c, d = (self.xmax, self.ymin)
# a, b = (self.xmin, self.ymax)
# e, f = (geodict.xmin, geodict.ymax)
# if geodict.xmin == 180 and geodict.xmax < 0:
# e, f = (-geodict.xmin, geodict.ymax)
# else:
# e, f = (geodict.xmin, geodict.ymax)
# g, h = (geodict.xmax, geodict.ymin)
# inside_x = (e >= a and e < c) or (a >= e and a < c)
# inside_y = (h >= d and h < b) or (d >= h and d < f)
# # inside_x = geodict.xmin >= self._xmin and geodict.xmax <= self._xmax
# # inside_y = geodict.ymin >= self._ymin and geodict.ymax <= self._ymax
# if inside_x and inside_y:
# return True
# return False
def _inside_x(self, geodict):
inside_x = False
if np.abs((self.xmax + self.dx) - (self.xmin + 360)) < self.dx * 0.01:
return True
if not inside_x:
if (
self.xmin > geodict.xmin
and self.xmax > geodict.xmin
and self.xmax < self.xmin
):
inside_x = (self.xmin - 360) < geodict.xmin and self.xmax > geodict.xmax
if inside_x:
return True
if geodict.xmin < self.xmin and (geodict.xmin + 360) < self.xmax:
inside_x = True
else:
if self.xmax > self.xmin:
xmin_inside = (geodict.xmin + self.EPS) >= self._xmin
xmax_inside = (geodict.xmax - self.EPS) <= self._xmax
else:
xmin_inside = (geodict.xmin + self.EPS) >= self._xmin
xmax_inside = (geodict.xmax - self.EPS) <= self._xmax + 360
inside_x = xmin_inside and xmax_inside
return inside_x
def contains(self, geodict):
"""Determine if input geodict is completely outside this GeoDict.
:param geodict:
Input GeoDict object.
:returns:
True if input geodict is completely outside this GeoDict,
False if not.
"""
inside_x = self._inside_x(geodict)
ymin_inside = (geodict.ymin + self.EPS) >= self._ymin
ymax_inside = (geodict.ymax - self.EPS) <= self._ymax
inside_y = ymin_inside and ymax_inside
if inside_x and inside_y:
return True
return False
def asDict(self):
"""Return GeoDict object in dictionary representation.
:returns:
Dictionary containing the same fields as found in constructor.
"""
tdict = {}
tdict["xmin"] = self._xmin
tdict["xmax"] = self._xmax
tdict["ymin"] = self._ymin
tdict["ymax"] = self._ymax
tdict["dx"] = self._dx
tdict["dy"] = self._dy
tdict["ny"] = self._ny
tdict["nx"] = self._nx
return tdict
def __repr__(self):
"""Return a string representation of the object."""
rfmt = """Bounds: (%.4f,%.4f,%.4f,%.4f)\nDims: (%.4f,%.4f)\nShape: (%i,%i)"""
rtpl = (
self._xmin,
self._xmax,
self._ymin,
self._ymax,
self._dx,
self._dy,
self._ny,
self._nx,
)
return rfmt % rtpl
def copy(self):
"""Return an object that is a complete copy of this GeoDict.
:returns:
A GeoDict object whose elements (xmin, xmax, ...) are an exact copy of
this object's elements.
"""
if not hasattr(self.__class__, "projection"):
self._projection = self.DEFAULT_PROJ4
geodict = {
"xmin": self._xmin,
"xmax": self._xmax,
"ymin": self._ymin,
"ymax": self._ymax,
"dx": self._dx,
"dy": self._dy,
"ny": self._ny,
"nx": self._nx,
"projection": self._projection,
}
return GeoDict(geodict)
def __eq__(self, other):
"""Check for equality between one GeoDict object and another.
:param other:
Another GeoDict object.
:returns:
True when all GeoDict parameters are no more different than 1e-12, False
otherwise.
"""
if np.abs(self._xmin - other._xmin) > self.EPS:
return False
if np.abs(self._ymin - other.ymin) > self.EPS:
return False
if np.abs(self._xmax - other.xmax) > self.EPS:
return False
if np.abs(self._ymax - other.ymax) > self.EPS:
return False
if np.abs(self._dx - other.dx) > self.EPS:
return False
if np.abs(self._dy - other.dy) > self.EPS:
return False
if np.abs(self._ny - other.ny) > self.EPS:
return False
if np.abs(self._nx - other.nx) > self.EPS:
return False
return True
def getLatLon(self, row, col):
"""Return geographic coordinates (lat/lon decimal degrees) for given data row
and column.
:param row:
Row dimension index into internal data array.
:param col:
Column dimension index into internal data array.
:returns:
Tuple of latitude and longitude.
"""
scalar_types = (int, float)
if type(row) != type(col):
raise DataSetException("Input row/col types must be the same")
if not isinstance(row, scalar_types): # is this a sequence type thing?
if isinstance(row, np.ndarray):
if row.shape == () or col.shape == ():
raise DataSetException(
"Input row/col values must be scalars or dimensioned numpy "
"arrays."
)
else:
raise DataSetException(
"Input row/col values must be scalars or dimensioned numpy arrays."
)
ulx = self._xmin
uly = self._ymax
dx = self._dx
dy = self._dy
lon = ulx + col * dx
lat = uly - row * dy
return (lat, lon)
def getRowCol(self, lat, lon, returnFloat=False, intMethod="round"):
"""Return data row and column from given geographic coordinates (lat/lon decimal
degrees).
:param lat:
Input latitude.
:param lon:
Input longitude.
:param returnFloat:
Boolean indicating whether floating point row/col coordinates should be
returned.
:param intMethod:
String indicating the desired method by which floating point row/col values
should
be converted to integers. Choices are 'round' (default), 'floor', or 'ceil'.
:returns:
Tuple of row and column.
"""
scalar_types = (int, float)
if type(lat) != type(lon):
raise DataSetException("Input lat/lon types must be the same")
if not isinstance(lat, scalar_types): # is this a sequence type thing?
if isinstance(lat, np.ndarray):
if lat.shape == () or lon.shape == ():
raise DataSetException(
"Input lat/lon values must be scalars or dimensioned numpy "
"arrays."
)
else:
raise DataSetException(
"Input lat/lon values must be scalars or dimensioned numpy arrays."
)
if intMethod not in ["round", "floor", "ceil"]:
raise DataSetException("intMethod %s is not supported." % intMethod)
ulx = self._xmin
uly = self._ymax
dx = self._dx
dy = self._dy
# check to see if we're in a scenario where the grid crosses the meridian
if self._xmax < ulx and np.any(lon < 0):
if not isinstance(lat, scalar_types):
lon[lon < 0] += 360
else:
lon += 360
col = (lon - ulx) / dx
row = (uly - lat) / dy
if returnFloat:
return (row, col)
if intMethod == "round":
return (np.round(row).astype(int), np.round(col).astype(int))
elif intMethod == "floor":
return (np.floor(row).astype(int), np.floor(col).astype(int))
else:
return (np.ceil(row).astype(int), np.ceil(col).astype(int))
# define setter and getter methods for all of the geodict parameters
@property
def xmin(self):
"""Get xmin value.
:returns:
xmin value.
"""
return self._xmin
@property
def xmax(self):
"""Get xmin value.
:returns:
xmin value.
"""
return self._xmax
@property
def ymin(self):
"""Get xmax value.
:returns:
xmax value.
"""
return self._ymin
@property
def ymax(self):
"""Get xmax value.
:returns:
xmax value.
"""
return self._ymax
@property
def dx(self):
"""Get dx value.
:returns:
dx value.
"""
return self._dx
@property
def dy(self):
"""Get dy value.
:returns:
dy value.
"""
return self._dy
@property
def ny(self):
"""Get ny value.
:returns:
ny value.
"""
return self._ny
@property
def nx(self):
"""Get nx value.
:returns:
nx value.
"""
return self._nx
@property
def projection(self):
"""Get projection Proj4 string.
:returns:
Valid Proj4 string.
"""
return self._projection
@property
def nodata(self):
"""Get projection Proj4 string.
:returns:
Valid nodata value.
"""
return self._nodata
@xmin.setter
def xmin(self, value):
"""Set xmin value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._xmin = value
self.validate()
@xmax.setter
def xmax(self, value):
"""Set xmax value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._xmax = value
self.validate()
@ymin.setter
def ymin(self, value):
"""Set ymin value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._ymin = value
self.validate()
@ymax.setter
def ymax(self, value):
"""Set ymax value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._ymax = value
self.validate()
@dx.setter
def dx(self, value):
"""Set dx value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._dx = value
self.validate()
@dy.setter
def dy(self, value):
"""Set dy value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._dy = value
self.validate()
@ny.setter
def ny(self, value):
"""Set ny value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._ny = value
self.validate()
@nx.setter
def nx(self, value):
"""Set nx value, re-validate object.
:param value:
Value to set.
:raises DataSetException:
When validation fails.
"""
self._nx = value
self.validate()
@nodata.setter
def nodata(self, value):
"""Set nodata value.
:param value:
Value to set.
"""
self._nodata = value
def getDeltas(self):
# handle cases where we're crossing the meridian from the eastern hemisphere
# to the western
if self._xmin > self._xmax:
txmax = self._xmax + 360.0
else:
txmax = self._xmax
# try calculating xmax from xmin, dx, and nx
xmax = self._xmin + self._dx * (self._nx - 1)
# dxmax = np.abs(xmax - txmax)
dxmax = np.abs(float(xmax) / txmax - 1.0)
# try calculating dx from bounds and nx
dx = np.abs((txmax - self._xmin) / (self._nx - 1))
# ddx = np.abs((dx - self._dx))
ddx = np.abs(float(dx) / self._dx - 1.0)
# try calculating ymax from ymin, dy, and ny
ymax = self._ymin + self._dy * (self._ny - 1)
# dymax = np.abs(ymax - self._ymax)
dymax = np.abs(float(ymax) / self._ymax - 1.0)
# try calculating dx from bounds and nx
dy = np.abs((self._ymax - self._ymin) / (self._ny - 1))
# ddy = np.abs(dy - self._dy)
ddy = np.abs(float(dy) / self._dy - 1.0)
return (dxmax, ddx, dymax, ddy)
def validate(self, adjust):
# dxmax,ddx,dymax,ddy = self.getDeltas()
if adjust == "bounds":
self._xmax = self._xmin + self._dx * (self._nx - 1)
self._ymin = self._ymax - self._dy * (self._ny - 1)
elif adjust == "res":
self._dx = (self._xmax - self._xmin) / (self._nx - 1)
if self._dx < 0:
if self._xmin > self._xmax:
self._dx = ((self._xmax + 360) - self._xmin) / (self._nx - 1)
self._dy = (self._ymax - self._ymin) / (self._ny - 1)
else:
raise DataSetException('Unsupported adjust option "%s"' % adjust)
if (self._xmax - self._xmin) < 359: # TODO - think about this more!
if self._xmax > 180:
self._xmax -= 360
if self._xmin < -180:
self._xmin += 360
# if self._xmin == 180.0 and self._xmax < 0:
# self._xmin = -180.0
def get_affine(src):
if rasterio.__version__ < "1.0.0":
affine = src.affine
else:
affine = src.transform
return affine
def geodict_from_src(src):
"""Return a geodict from a dataset object.
Args:
src (DatasetReader): Open rasterio DatasetReader object.