1 #!/usr/bin/env python3

2 

3 import argparse

4 import numpy as np

5 import sys

6 

7 def centroidnp(points):

8  lon_coords = [p[0] for p in points]

9  lat_coords = [p[1] for p in points]

10  length = len(points)

11  centroid_lon = sum(lon_coords)/length

12  centroid_lat = sum(lat_coords)/length

13  return [centroid_lon, centroid_lat]

14 

15 def monotonic(x):

16  dx = np.diff(x)

17  return np.all(dx <= 0) or np.all(dx >= 0)

18 

19 __version__ = "1.0-20211104"

20 

21 parser = argparse.ArgumentParser(prog='sort_gring',formatter_class=argparse.RawTextHelpFormatter,

22  description='This script sorts GRing points to ensure counter-clockwise order',epilog='''

23 EXIT Status:

24 0 : All is well in the world

25 1 : Dunno, something horrible occurred

26 110 : No reordering necessary''')

27 parser.add_argument('--version', action='version', version='%(prog)s ' + __version__)

28 parser.add_argument('-i','--ifile', type=str, metavar="FILE",default=None, help="name of file with gring data from l1info")

29 parser.add_argument('--lon', type=str, default=None, help='list of gring longitudes')

30 parser.add_argument('--lat', type=str, default=None, help="list of gring latitudes")

31 parser.add_argument('--verbose', '-v', action='store_true', default=False)

32 

33 args = parser.parse_args()

34 

35 gringpointlongitude = None

36 gringpointlatitude = None

37 

38 if args.ifile:

39  with open(args.ifile) as ifile:

40  lines = ifile.readlines()

41  lines = filter(str.rstrip, lines)

42  for line in lines:

43  if "=" in line:

44  key, value = line.split('=')

45  values = np.array(value.split(',')).astype(float)

46  if 'longitude' in key:

47  gringpointlongitude = np.squeeze(np.radians(np.array([values])))

48  if 'latitude' in key:

49  gringpointlatitude = np.squeeze(np.radians(np.array([values])))

50 

51 if args.lon:

52  value = args.lon

53  values = np.array(value.split(',')).astype(float)

54  gringpointlongitude = np.squeeze(np.radians(np.array([values])))

55 

56 if args.lat:

57  value = args.lat

58  values = np.array(value.split(',')).astype(float)

59  gringpointlatitude = np.squeeze(np.radians(np.array([values])))

60 

61 if (gringpointlongitude is not None) and (gringpointlatitude is not None):

62  if gringpointlongitude.shape != gringpointlatitude.shape:

63  sys.exit("lat and lon elements of different size")

64 

65  # make a list of tuples for the points

66  points = list(zip(gringpointlongitude, gringpointlatitude))

67  # find the middle of the polygon

68  center = centroidnp(points)

69  # get the distance from each point to the center

70  distance = np.atleast_1d(points) - np.atleast_1d(center)

71  # compute an angle from these

72  angles = np.arctan2(distance[:,0],distance[:,1])

73  # convert the angles into degrees in the range of 0-360

74  deg_angles = np.degrees(angles)

75  idx = np.where(deg_angles < 0)

76  deg_angles[idx] += 360

77  # sort the angles and reverse since we want counter-clockwise

78  gringsequence = np.argsort(deg_angles, axis=0)[::-1]

79 

80  if monotonic(gringsequence):

81  if args.verbose:

82  print("already ordered appropriately")

83  sys.exit(110)

84  else:

85  sorted_lon = np.array2string(np.take_along_axis(np.degrees(gringpointlongitude), gringsequence, axis=0), precision=5, separator=',')

86  sorted_lat = np.array2string(np.take_along_axis(np.degrees(gringpointlatitude), gringsequence, axis=0), precision=5, separator=',')

87  if args.verbose:

88  print("new world order: {}".format(gringsequence))

89 

90  print("gringpointlongitude={}".format(sorted_lon[1:-1]))

91  print("gringpointlatitude={}".format(sorted_lat[1:-1]))

92 else:

93  parser.print_help()