Skip to content

Commit

Permalink
Fix last minor stuffs + Code linting + Try to make tests run faster +…
Browse files Browse the repository at this point in the history 
… Bump patch nb version for release on pypi
  • Loading branch information
mz committed Sep 15, 2016
1 parent 442cbf3 commit c2f757e
Show file tree
Hide file tree
Showing 3 changed files with 78 additions and 61 deletions.
2 changes: 1 addition & 1 deletion smoomapy/__init__.py
View file
Expand Up @@ -4,6 +4,6 @@
"""
Smoomapy : make smoothed map in a python environnement.
"""
__version__ = "0.1.1"
__version__ = "0.1.2"

from .core import *
84 changes: 40 additions & 44 deletions smoomapy/core.py
View file
@@ -1,16 +1,14 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
More or less a python port of Stewart method
from R SpatialPositon package (https://github.com/Groupe-ElementR/SpatialPosition/)
Allow to set a desired number of class and choose discretization method
or directly set some custom breaks values.
More or less a python port of Stewart method from R SpatialPositon package
(https://github.com/Groupe-ElementR/SpatialPosition/).
@author: mthh
"""
import numpy as np
import pyproj
from scipy.interpolate import griddata as scipy_griddata, Rbf
from scipy.interpolate import griddata as scipy_griddata
from matplotlib.mlab import griddata as mlab_griddata
from matplotlib.pyplot import contourf
from shapely.geometry import Polygon, MultiPolygon, Point
Expand Down Expand Up @@ -47,16 +45,17 @@ def quick_stewart(input_geojson_points, variable_name, span,
beta : float
The beta!
typefct : str, optionnal
The type of function in {"exponential", "pareto"} (default: "exponential")
The type of function in {"exponential", "pareto"} (default: "exponential").
nb_class : int, default None
The number of class, if unset will most likely be 8
(default: None)
resolution : int, optionnal
The resolution to use (in unit of the input file), if not set a resolution
will be used in order to make a grid containing around 7560 pts
(default: None).
The resolution to use (in unit of the input file), if not set a default
resolution will be used in order to make a grid containing around
7500 pts (default: None).
mask : str, optionnal
Path to the file (Polygons only) to use as clipping mask (default: None).
Path to the file (Polygons only) to use as clipping mask,
can also be a GeoDataFrame (default: None).
user_defined_breaks : list or tuple, optionnal
A list of ordered break to use to construct the contours
(override `nb_class` value if any, default: None).
Expand All @@ -65,8 +64,8 @@ def quick_stewart(input_geojson_points, variable_name, span,
computed from this variable will be will be used as to divide
values computed from the first variable (default: None)
output : string, optionnal
The type of output expected (not case-sensitive) in {"GeoJSON", "GeoDataFrame"}
(default: "GeoJSON")
The type of output expected (not case-sensitive)
in {"GeoJSON", "GeoDataFrame"} (default: "GeoJSON").
Returns
-------
Expand Down Expand Up @@ -225,14 +224,8 @@ def hav_dist(locs1, locs2, k=np.pi/180):
cos_lat2 = np.cos(locs2[..., 0])
cos_lat_d = np.cos(locs1[..., 0] - locs2[..., 0])
cos_lon_d = np.cos(locs1[..., 1] - locs2[..., 1])
return 6367000 * np.arccos(cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))

#def check_bounds(xmin, ymin, xmax, ymax):
# if ymin < -9072187:
# ymin = -9072187
# if ymax > 9072187.8573143743:
# ymax = 9072187.8573143743
# return xmin, ymin, xmax, ymax
return 6367000 * np.arccos(
cos_lat_d - cos_lat1 * cos_lat2 * (1 - cos_lon_d))


def isopoly_to_gdf(collec_poly, levels, field_name="levels"):
Expand All @@ -248,7 +241,8 @@ def isopoly_to_gdf(collec_poly, levels, field_name="levels"):
levels : array-like
The value to use as attributes for the constructed GeoDataFrame.
field_name : str
The name of the field to be fill by `levels` variable (default: "levels")
The name of the field to be fill by values contained in
`levels` variable (default: "levels").
Returns
-------
Expand Down Expand Up @@ -276,12 +270,9 @@ def isopoly_to_gdf(collec_poly, levels, field_name="levels"):
polygons.append(mpoly[0])
data.append(levels[i])

# if len(data) == len(polygons):
return GeoDataFrame(geometry=polygons,
data=data,
columns=[field_name])
# else:
# return GeoDataFrame(geometry=polygons)


class Idx:
Expand Down Expand Up @@ -311,11 +302,11 @@ class SmoothStewart:
beta : float
The beta!
typefct : str, optionnal
The type of function in {"exponential", "pareto"} (default: "exponential")
The type of function in {"exponential", "pareto"} (default: "exponential").
resolution : int, optionnal
The resolution to use (in unit of the input file), if not set a resolution
will be used in order to make a grid containing around 7560 pts
(default: None).
The resolution to use (in unit of the input file), if not set a default
resolution will be used in order to make a grid containing around
7500 pts (default: None).
mask : str, optionnal
Path to the file (Polygons only) to use as clipping mask (default: None).
variable_name2 : str, optionnal
Expand Down Expand Up @@ -346,7 +337,9 @@ def __init__(self, input_layer, variable_name, span, beta,
self.longlat = kwargs.get("distGeo", kwargs.get("longlat", True))
self.gdf = input_layer if isinstance(input_layer, GeoDataFrame) \
else GeoDataFrame.from_file(input_layer)
self.proj_robinson = """+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"""
self.proj_robinson = (
"""+proj=robin +lon_0=0 +x_0=0 +y_0=0 """
"""+ellps=WGS84 +datum=WGS84 +units=m +no_defs""")
self.proj_nat_earth = """+proj=natearth"""
self.info = (
'SmoothStewart - variable : {}{} ({} features)\n'
Expand All @@ -365,13 +358,7 @@ def __init__(self, input_layer, variable_name, span, beta,
self.mask = GeoDataFrame.from_file(mask)

self.mask.to_crs(crs="+proj=natearth", inplace=True)

if len(set(self.mask.type)
.intersection({"Polygon", "MultiPolygon"})) > 0 \
and self.gdf.crs == self.mask.crs:
self.use_mask = True
else:
self.use_mask = False
self.check_mask()
else:
self.use_mask = False

Expand Down Expand Up @@ -437,7 +424,6 @@ def compute_pot(self, variable_name, span, beta,
[(knownpts.total_bounds[2] - knownpts.total_bounds[0])/10,
(knownpts.total_bounds[3] - knownpts.total_bounds[1])/10])
tmp = span * 1.5 if tmp < span * 1.5 else tmp
# bounds = check_bounds(*knownpts.buffer(tmp).total_bounds)
bounds = knownpts.buffer(tmp).total_bounds

self.unknownpts, self.shape = make_regular_points(bounds, resolution) \
Expand Down Expand Up @@ -545,6 +531,14 @@ def define_levels(self, nb_class, disc_func):

return levels

def check_mask(self):
if len(set(self.mask.type)
.intersection({"Polygon", "MultiPolygon"})) > 0 \
and self.gdf.crs == self.mask.crs:
self.use_mask = True
else:
self.use_mask = False

def render(self, nb_class=8, disc_func=None, user_defined_breaks=None,
func_grid="scipy", output="GeoJSON", new_mask=False):
"""
Expand All @@ -566,8 +560,12 @@ def render(self, nb_class=8, disc_func=None, user_defined_breaks=None,
"matplotlib" for matplotlib.mlab.griddata,
"rbf" for scipy.interpolate.rbf, default: "scipy").
output : string, optionnal
The type of output expected (not case-sensitive) in {"GeoJSON", "GeoDataFrame"}
(default: "GeoJSON").
The type of output expected (not case-sensitive)
in {"GeoJSON", "GeoDataFrame"} (default: "GeoJSON").
new_mask : str, optionnal
Use a new mask by giving the path to the file (Polygons only)
to use as clipping mask, can also be directly a GeoDataFrame
(default: False).
Returns
-------
Expand All @@ -586,10 +584,12 @@ def render(self, nb_class=8, disc_func=None, user_defined_breaks=None,
elif isinstance(new_mask, GeoDataFrame):
self.use_mask = True
self.mask = new_mask.to_crs(crs="+proj=natearth")
self.check_mask()
elif isinstance(new_mask, (str, BytesIO, StringIO)):
self.use_mask = True
self.mask = GeoDataFrame.from_file(
new_mask).to_crs(crs="+proj=natearth")
self.check_mask()

if func_grid == "scipy":
self.zi = scipy_griddata((self.x, self.y), pot,
Expand All @@ -607,10 +607,6 @@ def render(self, nb_class=8, disc_func=None, user_defined_breaks=None,
self.zi = mlab_griddata(self.x, self.y, pot,
self.xi, self.yi, interp='linear'
).round(8)
# elif func_grid == "rbf":
# rbf = Rbf(self.x, self.y, pot, epsilon=2)
# XI, YI = np.meshgrid(self.xi, self.yi)
# self.zi = rbf(XI, YI).round(8)
else:
raise ValueError("Invalid interpolation function name provided")

Expand Down Expand Up @@ -639,7 +635,7 @@ def render(self, nb_class=8, disc_func=None, user_defined_breaks=None,
res.loc[0:ix_max_ft, "geometry"] = res.geometry.buffer(
0).intersection(self.poly_max_extend.buffer(-0.1))
# Repair geometries if necessary :
if not all(t == "MultiPolygon" or t == "Polygon" for t in res.geom_type):
if not all(t in ("MultiPolygon", "Polygon") for t in res.geom_type):
res.loc[0:ix_max_ft, "geometry"] = \
[geom if geom.type in ("Polygon", "MultiPolygon")
else MultiPolygon(
Expand Down
53 changes: 37 additions & 16 deletions tests.py
View file
Expand Up @@ -18,15 +18,15 @@ def test_one_shot_stewart(self):
# Exports correctly to `bytes`:
res = quick_stewart(
"misc/nuts3_data.geojson", "pop2008",
span=65000, beta=2, resolution=55000, nb_class=8,
span=65000, beta=2, resolution=60000, nb_class=8,
mask="misc/nuts3_data.geojson")
self.assertIsInstance(res, bytes)

# Exports correctly to `GeoDataFrame`
# and respects the choosen number of class:
res = quick_stewart(
"misc/nuts3_data.geojson", "pop2008",
span=65000, beta=2, resolution=55000, nb_class=8,
span=65000, beta=2, resolution=60000, nb_class=8,
mask="misc/nuts3_data.geojson", output="GeoDataFrame")
self.assertIsInstance(res, GeoDataFrame)
self.assertEqual(len(res), 8)
Expand All @@ -48,7 +48,7 @@ def test_one_shot_stewart(self):
"pop2008",
span=65000,
beta=2,
resolution=55000,
resolution=60000,
user_defined_breaks=my_breaks,
mask="misc/nuts3_data.geojson",
output="GeoDataFrame")
Expand All @@ -61,8 +61,9 @@ def test_one_shot_stewart(self):
def test_object_stewart(self):
# Test the OO approach for building smoothed map with stewart potentials
StePot = SmoothStewart("misc/nuts3_data.geojson", "pop2008",
span=65000, beta=2, resolution=48000,
span=65000, beta=2, resolution=60000,
mask="misc/nuts3_data.geojson")

# Test using percentiles :
result = StePot.render(nb_class=10,
disc_func="percentiles",
Expand All @@ -81,26 +82,31 @@ def test_object_stewart(self):
# Test using somes already choosed break values :
my_breaks = [0, 1697631, 3395263, 5092894, 6790526,
8488157, 10185789, 11883420, 13581052]
result = StePot.render(nb_class=48, # bogus values as `nb_class` and
disc_func="foobar", # ... disc_func should be overrided
user_defined_breaks=my_breaks, # ... by the `user_defined_breaks` params
output="geodataframe") # ... and this is what we are testing here
result = StePot.render(
nb_class=48, # bogus values as `nb_class` and
disc_func="foobar", # ... disc_func should be overrided
user_defined_breaks=my_breaks, # ... by the `user_defined_breaks` params
output="geodataframe") # ... and this is what we are testing here

self.assertIsInstance(result, GeoDataFrame)
self.assertEqual(len(result), 8)
# Assert these break values were actually used :
for wanted_break, obtained_break in zip(my_breaks[1:-1], result["max"][:-1]):
self.assertAlmostEqual(wanted_break, obtained_break)

# Some tests on an other variables for testing another discretization method :
StePot = SmoothStewart("misc/nuts3_data.geojson", "gdppps2008",
span=65000, beta=2, resolution=48000,
mask="misc/nuts3_data.geojson")
# Using "head tail breaks" (should define automatically the number of class)
# Test again using another discretization method : "head tail breaks"
# (should define automatically the number of class)
result = StePot.render(nb_class=None,
disc_func="head_tail",
output="geodataframe")
self.assertIsInstance(result, GeoDataFrame)

# Test that the object has a nice representation :
a = str(StePot)
b = repr(StePot)
self.assertEqual(a, b)
self.assertIn("SmoothStewart - variable :", a)
self.assertIn("{} features".format(len(StePot.gdf)), a)

def test_object_stewart_two_var(self):
# Test the OO approach with two variables :
Expand Down Expand Up @@ -145,6 +151,7 @@ def test_from_gdf_with_new_mask(self):
result = StePot.render(5, output="Geodataframe",
new_mask="misc/nuts3_data.geojson")
self.assertIsInstance(result, GeoDataFrame)
self.assertEqual(StePot.use_mask, True)
self.assertEqual(len(result), 5)

# Or from a GeoDataFrame :
Expand All @@ -153,12 +160,24 @@ def test_from_gdf_with_new_mask(self):
result = StePot.render(6, output="Geodataframe",
new_mask=gdf)
self.assertIsInstance(result, GeoDataFrame)
self.assertEqual(StePot.use_mask, True)
self.assertEqual(len(result), 6)

# Nope, no mask :
result = StePot.render(5, output="Geodataframe",
new_mask=None)
self.assertIsInstance(result, GeoDataFrame)
self.assertEqual(StePot.use_mask, False)
self.assertEqual(len(result), 5)

# Test that it skips the mask parameter if the layer provided as a mask
# is not a Polygon/MultiPolygon layer :
gdf_mask = gdf[1:50].copy()
gdf_mask.geometry = gdf_mask.geometry.centroid
result = StePot.render(5, output="Geodataframe",
new_mask=gdf_mask)
self.assertIsInstance(result, GeoDataFrame)
self.assertEqual(StePot.use_mask, False)
self.assertEqual(len(result), 5)

def test_input_with_missing_values(self):
Expand All @@ -177,7 +196,7 @@ def test_wrong_dtype_missing_values(self):
gdf.loc[25:35, "pop2008"] = np.NaN
gdf.loc[0:len(gdf)-1, "pop2008"] = gdf["pop2008"].astype(str)
StePot = SmoothStewart(gdf, "gdppps2008",
span=65000, beta=2, resolution=48000,
span=65000, beta=2, resolution=60000,
mask="misc/nuts3_data.geojson")
result = StePot.render(9, "equal_interval", output="Geodataframe")
self.assertIsInstance(result, GeoDataFrame)
Expand All @@ -191,7 +210,8 @@ def test_errors(self):
typefct="abcdefg")

StePot = SmoothStewart("misc/nuts3_data.geojson", "gdppps2008",
span=65000, beta=2, resolution=48000)
span=65000, beta=2, resolution=60000)

# Test with a wrong discretization function name :
with self.assertRaises(ValueError):
StePot.render(9, "foo", output="Geodataframe")
Expand All @@ -210,8 +230,9 @@ def test_errors(self):

def test_mod_shape_interpolation_grid(self):
StePot = SmoothStewart("misc/nuts3_data.geojson", "pop2008",
span=80000, beta=2, resolution=75000,
span=65000, beta=2, resolution=75000,
mask="misc/nuts3_data.geojson")

# First rendering :
StePot.render(nb_class=8,
disc_func="percentiles",
Expand Down

0 comments on commit c2f757e

Please sign in to comment.