1""" 

2@file 

3@brief Geographic datasets. 

4""" 

5import math 

6import os 

7import pandas 

8from pyquickhelper.filehelper import un7zip_files 

9from .http_retrieve import download_data 

10 

11 

12def lambert932WGPS(lambertE, lambertN): 

13 """ 

14 Converts Lambert coordinates into longitude, latitude. 

15 Does not use :epkg:`pyproj`. 

16 

17 @param lambertE lambertE 

18 @param lambertN lambertN 

19 @return longitude, latitude 

20 """ 

21 

22 class constantes: 

23 GRS80E = 0.081819191042816 

24 LONG_0 = 3 

25 XS = 700000 

26 YS = 12655612.0499 

27 n = 0.7256077650532670 

28 C = 11754255.4261 

29 

30 delX = lambertE - constantes.XS 

31 delY = lambertN - constantes.YS 

32 gamma = math.atan(-delX / delY) 

33 R = math.sqrt(delX * delX + delY * delY) 

34 latiso = math.log(constantes.C / R) / constantes.n 

35 sinPhiit0 = math.tanh(latiso + constantes.GRS80E * 

36 math.atanh(constantes.GRS80E * math.sin(1))) 

37 sinPhiit1 = math.tanh(latiso + constantes.GRS80E * 

38 math.atanh(constantes.GRS80E * sinPhiit0)) 

39 sinPhiit2 = math.tanh(latiso + constantes.GRS80E * 

40 math.atanh(constantes.GRS80E * sinPhiit1)) 

41 sinPhiit3 = math.tanh(latiso + constantes.GRS80E * 

42 math.atanh(constantes.GRS80E * sinPhiit2)) 

43 sinPhiit4 = math.tanh(latiso + constantes.GRS80E * 

44 math.atanh(constantes.GRS80E * sinPhiit3)) 

45 sinPhiit5 = math.tanh(latiso + constantes.GRS80E * 

46 math.atanh(constantes.GRS80E * sinPhiit4)) 

47 sinPhiit6 = math.tanh(latiso + constantes.GRS80E * 

48 math.atanh(constantes.GRS80E * sinPhiit5)) 

49 

50 longRad = math.asin(sinPhiit6) 

51 latRad = gamma / constantes.n + constantes.LONG_0 / 180 * math.pi 

52 

53 longitude = latRad / math.pi * 180 

54 latitude = longRad / math.pi * 180 

55 

56 return longitude, latitude 

57 

58 

59def load_french_departements(cache=None): 

60 """ 

61 Returns the definition of French departments. 

62 

63 @param cache cache folder 

64 @return French departments as a :epkg:`GeoDataFrame` 

65 

66 .. exref:: 

67 :title: Loads French departments polygons 

68 

69 Simple example to retrieve French departements. 

70 

71 .. runpython:: 

72 :showcode: 

73 

74 from pyensae.datasource import load_french_departements 

75 df = load_french_departements() 

76 print(df.head(2).T) 

77 """ 

78 if cache is None: 

79 cache = '.' 

80 # delayed import 

81 import shapefile 

82 import geopandas 

83 from shapely.geometry import Polygon 

84 from shapely.ops import unary_union 

85 

86 name = "GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01.7z" 

87 try: 

88 download_data(name, whereTo=cache, 

89 website="https://wxs-telechargement.ign.fr/oikr5jryiph0iwhw36053ptm/telechargement/inspire/" 

90 "GEOFLA_THEME-DEPARTEMENTS_2015_2$GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01/file/") 

91 except Exception as e: 

92 # au cas le site n'est pas accessible 

93 download_data(name, website="xd", whereTo=cache) 

94 

95 full_name = os.path.join(cache, name) 

96 un7zip_files(full_name, where_to=os.path.join(cache, "shapefiles")) 

97 departements = os.path.join(cache, 'shapefiles/GEOFLA_2-1_DEPARTEMENT_SHP_LAMB93_FXX_2015-12-01/GEOFLA/1_DONNEES_LIVRAISON_2015/' 

98 'GEOFLA_2-1_SHP_LAMB93_FR-ED152/DEPARTEMENT/DEPARTEMENT.shp') 

99 

100 r = shapefile.Reader(departements) 

101 shapes = r.shapes() 

102 records = r.records() 

103 

104 polys = [] 

105 datas = [] 

106 for i, (record, shape) in enumerate(zip(records, shapes)): 

107 # coordinates in Lambert 93 

108 geo_points = [lambert932WGPS(x, y) for x, y in shape.points] 

109 rec = {} 

110 for k in ['CODE_DEPT', 'CODE_REG', 'CODE_CHF', 'ID_GEOFLA', 'NOM_CHF', 

111 'NOM_DEPT', 'NOM_REG', 'X_CENTROID', 'X_CHF_LIEU', 

112 'Y_CENTROID', 'Y_CHF_LIEU']: 

113 rec[k] = getattr(record, k, None) 

114 if len(shape.parts) == 1: 

115 # only one polygon 

116 poly = Polygon(geo_points) 

117 else: 

118 # merge them 

119 ind = list(shape.parts) + [len(shape.points)] 

120 pols = [Polygon(geo_points[ind[i]:ind[i + 1]]) 

121 for i in range(0, len(shape.parts))] 

122 try: 

123 poly = unary_union(pols) 

124 except Exception as e: 

125 poly = Polygon(geo_points) 

126 polys.append(poly) 

127 datas.append(rec) 

128 

129 geom = geopandas.GeoDataFrame(geometry=polys) 

130 data = pandas.DataFrame(datas) 

131 return pandas.concat([geom, data], axis=1)