Coverage for src/pyensae/datasource/geodata.py: 96%
68 statements
« prev ^ index » next coverage.py v7.2.7, created at 2023-07-03 02:16 +0200
« prev ^ index » next coverage.py v7.2.7, created at 2023-07-03 02:16 +0200
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
12def lambert932WGPS(lambertE, lambertN):
13 """
14 Converts Lambert coordinates into longitude, latitude.
15 Does not use :epkg:`pyproj`.
17 @param lambertE lambertE
18 @param lambertN lambertN
19 @return longitude, latitude
20 """
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
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))
50 longRad = math.asin(sinPhiit6)
51 latRad = gamma / constantes.n + constantes.LONG_0 / 180 * math.pi
53 longitude = latRad / math.pi * 180
54 latitude = longRad / math.pi * 180
56 return longitude, latitude
59def load_french_departements(cache=None):
60 """
61 Returns the definition of French departments.
63 @param cache cache folder
64 @return French departments as a :epkg:`GeoDataFrame`
66 .. exref::
67 :title: Loads French departments polygons
69 Simple example to retrieve French departements.
71 .. runpython::
72 :showcode:
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
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)
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')
100 r = shapefile.Reader(departements)
101 shapes = r.shapes()
102 records = r.records()
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)
129 geom = geopandas.GeoDataFrame(geometry=polys)
130 data = pandas.DataFrame(datas)
131 return pandas.concat([geom, data], axis=1)