Coverage for hml/rasterizePolygon.py: 3%

1#!/usr/bin/env python3 

2 

3# Define function ... 

4def rasterizePolygon(poly, /, *, px = 1024.0): 

5 """ 

6 Rasterize a [Multi]Polygon. 

7 

8 Arguments: 

9 poly -- a shapely.geometry.[multi]polygon.[Multi]Polygon 

10 

11 Keyword arguments: 

12 px -- pixel size (default 1024.0) 

13 

14 Note: 

15 This function only works for [Multi]Polygons that solely exist in the 

16 (positive, positive) quadrant. 

17 """ 

18 

19 # Import standard modules ... 

20 import math 

21 

22 # Import special modules ... 

23 try: 

24 import numpy 

25 except: 

26 raise Exception("\"numpy\" is not installed; run \"pip install --user numpy\"") from None 

27 try: 

28 import shapely 

29 import shapely.geometry 

30 except: 

31 raise Exception("\"shapely\" is not installed; run \"pip install --user Shapely\"") from None 

32 

33 # Check argument ... 

34 if not isinstance(poly, shapely.geometry.polygon.Polygon): 

35 if not isinstance(poly, shapely.geometry.multipolygon.MultiPolygon): 

36 raise TypeError("\"poly\" is not a shapely.geometry.[multi]polygon.[Multi]Polygon") 

37 

38 # Find bounding pixel indices in the global grid ... 

39 ix1 = math.floor(poly.bounds[0] / px) 

40 iy1 = math.floor(poly.bounds[1] / px) 

41 ix2 = math.ceil(poly.bounds[2] / px) 

42 iy2 = math.ceil(poly.bounds[3] / px) 

43 

44 # Find extent of the local grid ... 

45 nx = ix2 - ix1 

46 ny = iy2 - iy1 

47 

48 # Initialize local grid ... 

49 localGrid = numpy.zeros((ny, nx), dtype = numpy.float32) # [m2] 

50 

51 # Loop over x-axis ... 

52 for ix in range(nx): 

53 # Create short-hands ... 

54 xmin = float(ix1 + ix) * px # [m] 

55 xmax = float(ix1 + ix + 1) * px # [m] 

56 

57 # Loop over y-axis ... 

58 for iy in range(ny): 

59 # Create short-hands ... 

60 ymin = float(iy1 + iy) * px # [m] 

61 ymax = float(iy1 + iy + 1) * px # [m] 

62 

63 # Create a counter-clockwise polygon of the pixel, find its 

64 # intersection with the [Multi]Polygon and add the area to the local 

65 # grid ... 

66 localGrid[iy, ix] += poly.intersection( 

67 shapely.geometry.polygon.Polygon( 

68 [ 

69 (xmin, ymin), 

70 (xmax, ymin), 

71 (xmax, ymax), 

72 (xmin, ymax), 

73 (xmin, ymin), 

74 ] 

75 ) 

76 ).area # [m2] 

77 

78 # Return answer ... 

79 return ix1, iy1, localGrid