From e5c9379262b76f46c4da2b8e5f984929859d245d Mon Sep 17 00:00:00 2001 From: Anand Patil Date: Thu, 3 Mar 2011 12:22:18 +0000 Subject: [PATCH] Rastervals_in_unit is now breaking down multipolygons. --- map_utils/shapefile_utils.py | 17 ++++++++++++----- 1 file changed, 12 insertions(+), 5 deletions(-) diff --git a/map_utils/shapefile_utils.py b/map_utils/shapefile_utils.py index f8a84f5..4741b37 100755 --- a/map_utils/shapefile_utils.py +++ b/map_utils/shapefile_utils.py @@ -169,6 +169,11 @@ def multipoly_sample(n, mp, test=None, verbose=0): return lons, lats def rastervals_in_unit(unit, lon_min, lat_min, cellsize, data, view='y-x+'): + + # For some reason it is MUCH faster to do it this way. + if isinstance(unit, geometry.multipolygon.MultiPolygon): + return np.hstack(*[rastervals_in_unit(g,lon_min,lat_min,cellsize,data,view) for g in unit.geoms]) + llc = unit.bounds[:2] urc = unit.bounds[2:] @@ -177,9 +182,9 @@ def rastervals_in_unit(unit, lon_min, lat_min, cellsize, data, view='y-x+'): lat_min_in = int((llc[1]-lat_min)/cellsize) lat_max_in = int((urc[1]-lat_min)/cellsize)+1 - data_boxed = grid_convert(data,view,'x+y+')[lon_min_in:lon_max_in+1,lat_min_in:lat_max_in+1] - flat_shape = (-1,)+data_boxed.shape[2:] - data_boxed = data_boxed.reshape(flat_shape) + data_boxed_ = grid_convert(data,view,'x+y+')[lon_min_in:lon_max_in+1,lat_min_in:lat_max_in+1] + flat_shape = (-1,)+data_boxed_.shape[2:] + data_boxed = data_boxed_.reshape(flat_shape) lon_boxed = np.arange(lon_min_in, lon_max_in+1)*cellsize + lon_min lat_boxed = np.arange(lat_min_in, lat_max_in+1)*cellsize + lat_min @@ -188,8 +193,10 @@ def rastervals_in_unit(unit, lon_min, lat_min, cellsize, data, view='y-x+'): mlat_boxed = grid_convert(mlat_boxed,'y+x+','x+y+') p=[geom.Point(*pair) for pair in zip(mlon_boxed.ravel(), mlat_boxed.ravel())] - out = data_boxed[np.where([unit.contains(p_) for p_ in p])] - + p_in = iterops.contains(unit, p, True) + i_in = [p.index(p_in_) for p_in_ in p] + out = data_boxed[i_in] + return out def unit_to_grid(unit, lon_min, lat_min, cellsize):