Skip to content

Commit

Permalink
various fixes for weight & voronoi
Browse files Browse the repository at this point in the history
  • Loading branch information
snowman2 committed Oct 7, 2016
1 parent 08289d1 commit 3944fcd
Show file tree
Hide file tree
Showing 2 changed files with 22 additions and 21 deletions.
17 changes: 6 additions & 11 deletions RAPIDpy/gis/voronoi.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,14 +22,10 @@ def _get_voronoi_centroid_array(lsm_lat_array, lsm_lon_array, extent):
This function generates a voronoi centroid point
list from arrays of latitude and longitude
"""
if extent:
YMin = extent[2]
YMax = extent[3]
XMin = extent[0]
XMax = extent[1]
else:
lsm_lat_list = lsm_lat_array
lsm_lon_list = lsm_lon_array
YMin = extent[2]
YMax = extent[3]
XMin = extent[0]
XMax = extent[1]

ptList = []
if (lsm_lat_array.ndim == 2) and (lsm_lon_array.ndim == 2):
Expand All @@ -51,7 +47,6 @@ def _get_voronoi_centroid_array(lsm_lat_array, lsm_lon_array, extent):

lsm_lat_list = lsm_lat_array[lsm_lat_indices,:][:,lsm_lon_indices]
lsm_lon_list = lsm_lon_array[lsm_lat_indices,:][:,lsm_lon_indices]

# Create a list of geographic coordinate pairs
for i in range(len(lsm_lat_indices)):
for j in range(len(lsm_lon_indices)):
Expand Down Expand Up @@ -174,8 +169,8 @@ def pointsToVoronoiGridShapefile(lat, lon, vor_shp_path, extent=None):
ring.AddPoint(loopLon,loopLat)
poly.AddGeometry(ring)
feat = ogr.Feature(layerDefn)
feat.SetField('GRID_LON', voronoi_centroid[0])
feat.SetField('GRID_LAT', voronoi_centroid[1])
feat.SetField('GRID_LON', float(voronoi_centroid[0]))
feat.SetField('GRID_LAT', float(voronoi_centroid[1]))
feat.SetGeometry(poly)
layer.CreateFeature(feat)
feat = poly = ring = None
Expand Down
26 changes: 16 additions & 10 deletions RAPIDpy/gis/weight.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,11 +53,11 @@ def _get_lat_lon_indices(lsm_lat_array, lsm_lon_array, lat, lon):
"""
if lsm_lat_array.ndim == 2 and lsm_lon_array.ndim == 2:
lsm_lat_indices_from_lat, lsm_lon_indices_from_lat = np.where((lsm_lat_array == lat))
lsm_lat_indices_from_lon, lsm_lon_indices_from_lon = np.where((lsm_lon_array >= lon))
lsm_lat_indices_from_lon, lsm_lon_indices_from_lon = np.where((lsm_lon_array == lon))

index_lsm_grid_lat = np.intersect1d(lsm_lat_indices_from_lat, lsm_lat_indices_from_lon)[0]
index_lsm_grid_lon = np.intersect1d(lsm_lon_indices_from_lat, lsm_lon_indices_from_lon)[0]

elif lsm_lat_array.ndim == 1 and lsm_lon_array.ndim == 1:
index_lsm_grid_lon = np.where(lsm_lon_array == lon)[0][0]
index_lsm_grid_lat = np.where(lsm_lat_array == lat)[0][0]
Expand All @@ -82,6 +82,11 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,
Create Weight Table for Land Surface Model Grids
"""
time_start_all = datetime.utcnow()

if lsm_grid_lat.ndim == 3 and lsm_grid_lon.ndim == 3:
#assume first dimension is time
lsm_grid_lat = lsm_grid_lat[0]
lsm_grid_lon = lsm_grid_lon[0]

print("Generating LSM Grid Thiessen Array ...")
if file_geodatabase:
Expand All @@ -106,10 +111,10 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,
lsm_grid_feature_list = pointsToVoronoiGridArray(lsm_grid_lat, lsm_grid_lon, extent)

##COMMENTED LINES FOR TESTING
##import os
##from .voronoi import pointsToVoronoiGridShapefile
##vor_shp_path = os.path.join(os.path.dirname(in_catchment_shapefile), "test_grid.shp")
##pointsToVoronoiGridShapefile(lsm_grid_lat, lsm_grid_lon, vor_shp_path, extent)
# import os
# from .voronoi import pointsToVoronoiGridShapefile
# vor_shp_path = os.path.join(os.path.dirname(in_catchment_shapefile), "test_grid.shp")
# pointsToVoronoiGridShapefile(lsm_grid_lat, lsm_grid_lon, vor_shp_path, extent)

time_end_lsm_grid_thiessen = datetime.utcnow()
print(time_end_lsm_grid_thiessen - time_start_all)
Expand Down Expand Up @@ -159,7 +164,6 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,

for rapid_connect_rivid in rapid_connect_rivid_list:
intersect_grid_info_list = []

try:
catchment_pos = np.where(catchment_rivid_list==rapid_connect_rivid)[0][0]
except IndexError:
Expand All @@ -172,7 +176,8 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,
#make sure coordinates are geographic
if proj_transform:
feat_geom.Transform(proj_transform)
catchment_polygon = shapely_loads(feat_geom.ExportToWkb())
catchment_polygon = shapely_loads(feat_geom.ExportToWkb())

for sub_lsm_grid_pos in rtree_idx.intersection(catchment_polygon.bounds):
if catchment_polygon.intersects(lsm_grid_feature_list[sub_lsm_grid_pos]['polygon']):
intersect_poly = catchment_polygon.intersection(lsm_grid_feature_list[sub_lsm_grid_pos]['polygon'])
Expand All @@ -181,7 +186,7 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,
poly_area = get_poly_area_geo(intersect_poly)
else:
poly_area = float(catchment_polygon.GetFeature(area_id))*intersect_poly.area/catchment_polygon.area

index_lsm_grid_lat, index_lsm_grid_lon = _get_lat_lon_indices(lsm_grid_lat, lsm_grid_lon,
lsm_grid_feature_list[sub_lsm_grid_pos]['lat'],
lsm_grid_feature_list[sub_lsm_grid_pos]['lon'])
Expand All @@ -191,6 +196,7 @@ def RTreeCreateWeightTable(lsm_grid_lat, lsm_grid_lon,
'lsm_grid_lon': lsm_grid_feature_list[sub_lsm_grid_pos]['lon'],
'index_lsm_grid_lon': index_lsm_grid_lon,
'index_lsm_grid_lat': index_lsm_grid_lat})

npoints = len(intersect_grid_info_list)
for intersect_grid_info in intersect_grid_info_list:
connectwriter.writerow([intersect_grid_info['rivid'],
Expand Down Expand Up @@ -272,7 +278,7 @@ def CreateWeightTableLDAS(in_ldas_nc,
file_geodatabase=None):

"""
Create Weight Table for NLDAS, GLDAS grids as well as for 2D WRF, Joules, or LIS Grids
Create Weight Table for NLDAS, GLDAS grids as well as for 2D Joules, or LIS Grids
Args:
in_ldas_nc(str): Path to the land surface model NetCDF grid.
Expand Down

0 comments on commit 3944fcd

Please sign in to comment.