johnziebro
commented
2 years ago In attempting to create a reproducable example I ran into another issue.
ValueError: Must have equal len keys and value when setting with an iterable
import requests
import geopandas as gpd
from tobler.dasymetric import extract_raster_features, masked_area_interpolate"""
# commented out as uneeded for example since data is loaded from gist
def miles_radius_by_coords(x_lon:float, y_lat:float, miles:float, gdf:gpd.geodataframe.GeoDataFrame, reset_index:bool=False) -> List[gpd.geodataframe.GeoDataFrame]:
#Provided a point in epsg:4326, buffers by miles and return intersecting polygons of provided gdf.
gdf_crs = gdf.crs
gdf_columns = gdf.columns.tolist()
window = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=[x_lon], y=[y_lat], crs="epsg:4326").to_crs("esri:102008").buffer(miles * 1609.344))
gdf = gdf.to_crs("esri:102008").sjoin(window, how="inner", predicate="intersects")[gdf_columns].to_crs(gdf_crs)
return gdf.reset_index(drop=True) if reset_index else gdf
# determine nlcd crs
mask = "./nlcd_2019_land_cover_l48_20210604.img"
nlcd_crs = extract_raster_features(target, mask, pixel_values=[21, 22, 23, 24], collapse_values=True).crs
# subset geodataframes by 2 mile radius from location, setting crs to match mask's
miles = 2
coords = (-84.558735, 39.230207)
source = miles_radius_by_coords(*coords, miles, score_geo, reset_index=True).to_crs(nlcd_crs)
target = miles_radius_by_coords(*coords, miles, gdf_select_zctas, reset_index=True).to_crs(nlcd_crs)"""# designate geotiff path
mask = "./nlcd_2019_land_cover_l48_20210604.img"
# determine nlcd crs
mask = "./nlcd_2019_land_cover_l48_20210604.img"
nlcd_crs = extract_raster_features(target, mask, pixel_values=[21, 22, 23, 24], collapse_values=True).crs
# load data from gist json features
source_url = "https://gist.githubusercontent.com/johnziebro/d1bc6c59f55a8ba5bf6ead4cd3d2d5ee/raw/3e8370228056f05bd5515ea93f1fdaf8af62f778/source"
target_url = "https://gist.githubusercontent.com/johnziebro/0531432faa14d51073674d471ef99d95/raw/a80a8b2a61fffa8a56f2d4fb1721e2278a3dc41c/target"
source = gpd.GeoDataFrame.from_features(requests.get(source_url).json(), crs=nlcd_crs)
target = gpd.GeoDataFrame.from_features(requests.get(target_url).json(), crs=nlcd_crs)
# sanity check
display(source, source.isna().sum(), source.crs, target, target.isna().sum(), target.crs)
source.plot()
target.plot()Results display as expected
# define intensive and extensive variables
# intensive: do not depend on sample size, ex: percentages
# extensive: depend on sample size: ex: population
intensive = [
"ALAND_SMI_S",
"POP_DENSITY_SMI_S",
"SCORE",
"SCORE_POP",
"SCORE_VEH",
"SCORE_HOME",
"SCORE_FIN",
]
extensive = ["ALAND", "AWATER", "ALAND_SMI", "POP_DENSITY_SMI"]
# interpolate
results = masked_area_interpolate(raster=mask,
source_df=source,
target_df=target,
intensive_variables=intensive,
extensive_variables=extensive)
results
knaaptime
Transferred from comments in https://github.com/pysal/tobler/issues/136 per @martinfleis request.
I am encountering an error using masked_area_interpolate() between ~700 US Census tracts to ~200 US Census ZCTAs with nlcd_2019_land_cover as the mask.
Verified features can be extracted successfully from the NLCD geotiff using extract_raster_features().
Example:
ensure crs matches NLCD, leveraging developed_area from previous example
nlcd_crs = developed_area.crs.to_string() source.to_crs(nlcd_crs, inplace=True) target.to_crs(nlcd_crs, inplace=True)
attempt possible solution supplied at https://github.com/pysal/tobler/issues/136#issuecomment-786244084
source.geometry = source.buffer(0) target.geometry = target.buffer(0)
define intensive and extensive variables
intensive: do not depend on sample size, ex: percentages
extensive: depend on sample size: ex: population
intensive = [ "ALAND_SMI_SCALED", "POP_DENSITY_SMI_SCALED", ] extensive = ["ALAND", "POP_DENSITY_SMI"]
masked interpolation
results = masked_area_interpolate(raster=mask, source_df=source, target_df=target, intensive_variables=intensive, extensive_variables=extensive)