In [4]:
'''
Access to improved water sources with children and elderly overlay
'''
import os
import rioxarray as rxr
import geopandas as gpd
import pandas as pd
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
import matplotlib.pyplot as plt
import contextily as ctx
import rasterio
from rasterio.plot import show
import rasterio.mask
import matplotlib.pyplot as plt
os.chdir(r'C:/Users/jtrum/world_bank/data/')
In [5]:
'''
aoi = area of interest
hdsl = high resolution population density of children under 5
ihme = mean percent of population with access to improved water source
'''
aoi = gpd.read_file('aoiLuanda.geojson')
hdsl = rxr.open_rasterio('ago_children_under_five_2020.tif')
ihme = rxr.open_rasterio('IHME_LMIC_WASH_2000_2017_W_IMP_PERCENT_MEAN_2017_Y2020M06D02.TIF')
aoi_bound = aoi.geometry
hdsl_clip = hdsl.rio.clip(aoi_bound)
ihme_clip = ihme.rio.clip(aoi_bound)
In [6]:
print(type(hdsl_clip))
print(type(ihme_clip))
print(type(aoi))
<class 'xarray.core.dataarray.DataArray'> <class 'xarray.core.dataarray.DataArray'> <class 'geopandas.geodataframe.GeoDataFrame'>
In [21]:
hdsl_clip2 = hdsl_clip.squeeze()
hdsl_df = hdsl_clip2.to_dataframe(name='percent')
hdsl_df = hdsl_df['percent'].dropna()
ihme_clip2 = ihme_clip.squeeze()
ihme_df = ihme_clip2.to_dataframe(name='percent')
ihme_df = ihme_df[ihme_df.percent != -999999.0]
In [22]:
hdsl_df = hdsl_df.reset_index()
hdsl_df = hdsl_df.rename(columns={'x':'lon', 'y':'lat'})
hdsl_gdf = gpd.GeoDataFrame(hdsl_df, geometry=gpd.points_from_xy(hdsl_df.lon, hdsl_df.lat))
In [23]:
ihme_df =ihme_df.reset_index()
In [24]:
ihme_df
Out[24]:
| y | x | band | spatial_ref | percent | |
|---|---|---|---|---|---|
| 0 | -8.687501 | 13.437516 | 1 | 0 | 99.396202 |
| 1 | -8.687501 | 13.479183 | 1 | 0 | 99.813766 |
| 2 | -8.687501 | 13.520849 | 1 | 0 | 96.345879 |
| 3 | -8.729168 | 13.437516 | 1 | 0 | 99.989731 |
| 4 | -8.729168 | 13.479183 | 1 | 0 | 99.981781 |
| ... | ... | ... | ... | ... | ... |
| 111 | -9.229168 | 13.187516 | 1 | 0 | 33.793404 |
| 112 | -9.270834 | 13.145849 | 1 | 0 | 36.304882 |
| 113 | -9.270834 | 13.187516 | 1 | 0 | 36.074020 |
| 114 | -9.312501 | 13.145849 | 1 | 0 | 37.536205 |
| 115 | -9.312501 | 13.187516 | 1 | 0 | 33.094032 |
116 rows × 5 columns
In [25]:
ihme_df = ihme_df.rename(columns={'x':'lon', 'y':'lat'})
In [26]:
ihme_gdf = gpd.GeoDataFrame(ihme_df, geometry=(ihme_df.lon, ihme_df.lat))
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\core\indexes\base.py:3653, in Index.get_loc(self, key) 3652 try: -> 3653 return self._engine.get_loc(casted_key) 3654 except KeyError as err: File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\_libs\index.pyx:147, in pandas._libs.index.IndexEngine.get_loc() File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\_libs\index.pyx:153, in pandas._libs.index.IndexEngine.get_loc() TypeError: '(0 13.437516 1 13.479183 2 13.520849 3 13.437516 4 13.479183 ... 111 13.187516 112 13.145849 113 13.187516 114 13.145849 115 13.187516 Name: lon, Length: 116, dtype: float64, 0 -8.687501 1 -8.687501 2 -8.687501 3 -8.729168 4 -8.729168 ... 111 -9.229168 112 -9.270834 113 -9.270834 114 -9.312501 115 -9.312501 Name: lat, Length: 116, dtype: float64)' is an invalid key During handling of the above exception, another exception occurred: InvalidIndexError Traceback (most recent call last) c:\Users\jtrum\Desktop\maps\map_11_water-sources_children-elderly.ipynb Cell 9 line 1 ----> <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#Y115sZmlsZQ%3D%3D?line=0'>1</a> ihme_gdf = gpd.GeoDataFrame(ihme_df, geometry=(ihme_df.lon, ihme_df.lat)) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\geopandas\geodataframe.py:191, in GeoDataFrame.__init__(self, data, geometry, crs, *args, **kwargs) 183 if ( 184 hasattr(geometry, "crs") 185 and geometry.crs 186 and crs 187 and not geometry.crs == crs 188 ): 189 raise ValueError(crs_mismatch_error) --> 191 self.set_geometry(geometry, inplace=True, crs=crs) 193 if geometry is None and crs: 194 raise ValueError( 195 "Assigning CRS to a GeoDataFrame without a geometry column is not " 196 "supported. Supply geometry using the 'geometry=' keyword argument, " 197 "or by providing a DataFrame with column name 'geometry'", 198 ) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\geopandas\geodataframe.py:320, in GeoDataFrame.set_geometry(self, col, drop, inplace, crs) 318 else: 319 try: --> 320 level = frame[col] 321 except KeyError: 322 raise ValueError("Unknown column %s" % col) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\geopandas\geodataframe.py:1475, in GeoDataFrame.__getitem__(self, key) 1469 def __getitem__(self, key): 1470 """ 1471 If the result is a column containing only 'geometry', return a 1472 GeoSeries. If it's a DataFrame with any columns of GeometryDtype, 1473 return a GeoDataFrame. 1474 """ -> 1475 result = super().__getitem__(key) 1476 # Custom logic to avoid waiting for pandas GH51895 1477 # result is not geometry dtype for multi-indexes 1478 if ( 1479 pd.api.types.is_scalar(key) 1480 and key == "" (...) 1483 and not is_geometry_type(result) 1484 ): File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\core\frame.py:3761, in DataFrame.__getitem__(self, key) 3759 if self.columns.nlevels > 1: 3760 return self._getitem_multilevel(key) -> 3761 indexer = self.columns.get_loc(key) 3762 if is_integer(indexer): 3763 indexer = [indexer] File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\core\indexes\base.py:3660, in Index.get_loc(self, key) 3655 raise KeyError(key) from err 3656 except TypeError: 3657 # If we have a listlike key, _check_indexing_error will raise 3658 # InvalidIndexError. Otherwise we fall through and re-raise 3659 # the TypeError. -> 3660 self._check_indexing_error(key) 3661 raise File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\pandas\core\indexes\base.py:5737, in Index._check_indexing_error(self, key) 5733 def _check_indexing_error(self, key): 5734 if not is_scalar(key): 5735 # if key is not a scalar, directly raise an error (the code below 5736 # would convert to numpy arrays and raise later any way) - GH29926 -> 5737 raise InvalidIndexError(key) InvalidIndexError: (0 13.437516 1 13.479183 2 13.520849 3 13.437516 4 13.479183 ... 111 13.187516 112 13.145849 113 13.187516 114 13.145849 115 13.187516 Name: lon, Length: 116, dtype: float64, 0 -8.687501 1 -8.687501 2 -8.687501 3 -8.729168 4 -8.729168 ... 111 -9.229168 112 -9.270834 113 -9.270834 114 -9.312501 115 -9.312501 Name: lat, Length: 116, dtype: float64)
In [12]:
ihme_gdf
Out[12]:
| lat | lon | band | spatial_ref | percent | geometry | |
|---|---|---|---|---|---|---|
| 0 | -8.687501 | 13.437516 | 1 | 0 | 99.396202 | POINT (13.43752 -8.68750) |
| 1 | -8.687501 | 13.479183 | 1 | 0 | 99.813766 | POINT (13.47918 -8.68750) |
| 2 | -8.687501 | 13.520849 | 1 | 0 | 96.345879 | POINT (13.52085 -8.68750) |
| 3 | -8.729168 | 13.437516 | 1 | 0 | 99.989731 | POINT (13.43752 -8.72917) |
| 4 | -8.729168 | 13.479183 | 1 | 0 | 99.981781 | POINT (13.47918 -8.72917) |
| ... | ... | ... | ... | ... | ... | ... |
| 111 | -9.229168 | 13.187516 | 1 | 0 | 33.793404 | POINT (13.18752 -9.22917) |
| 112 | -9.270834 | 13.145849 | 1 | 0 | 36.304882 | POINT (13.14585 -9.27083) |
| 113 | -9.270834 | 13.187516 | 1 | 0 | 36.074020 | POINT (13.18752 -9.27083) |
| 114 | -9.312501 | 13.145849 | 1 | 0 | 37.536205 | POINT (13.14585 -9.31250) |
| 115 | -9.312501 | 13.187516 | 1 | 0 | 33.094032 | POINT (13.18752 -9.31250) |
116 rows × 6 columns
In [13]:
# Reproject the AOI to match the CRS of the raster datasets
aoi = aoi.to_crs(hdsl.crs)
# Get the geometry of the AOI
aoi_geometry = aoi.geometry.values[0]
In [14]:
# Clip the raster datasets to the AOI
hdsl_clip, hdsl_transform = rasterio.mask.mask(hdsl, shapes=[aoi_geometry], crop=True)
ihme_clip, ihme_transform = rasterio.mask.mask(ihme, shapes=[aoi_geometry], crop=True)
In [16]:
hdsl_clip_dataset = rasterio.open(
'hdsl_clip.tif',
'w',
driver='GTiff',
height=hdsl_clip.shape[1],
width=hdsl_clip.shape[2],
count=1,
dtype=str(hdsl_clip.dtype),
crs=hdsl.crs,
transform=hdsl_transform,
)
hdsl_clip_dataset.write(hdsl_clip[0], 1)
hdsl_clip_dataset.close()
In [17]:
# Create a new GeoDataFrame for the clipped AOI
aoi_clipped = gpd.GeoDataFrame({'geometry': [aoi_geometry]}, crs=hdsl.crs)
In [18]:
# Manually set CRS for GeoDataFrame
aoi_clipped = aoi_clipped.to_crs(epsg=3857)
In [21]:
ax = aoi.plot(facecolor="none",
edgecolor="red",
linewidth=2
)
ctx.add_basemap(ax,
crs=aoi.crs.to_string(),
source=ctx.providers.CartoDB.Voyager
)
--------------------------------------------------------------------------- CPLE_BaseError Traceback (most recent call last) File rasterio\\crs.pyx:775, in rasterio.crs.CRS.from_user_input() File rasterio\\_err.pyx:209, in rasterio._err.exc_wrap_ogrerr() CPLE_BaseError: OGR Error code 6 During handling of the above exception, another exception occurred: CRSError Traceback (most recent call last) c:\Users\jtrum\Desktop\maps\map_11_water-sources_children-elderly.ipynb Cell 11 line 5 <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=0'>1</a> ax = aoi.plot(facecolor="none", <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=1'>2</a> edgecolor="red", <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=2'>3</a> linewidth=2 <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=3'>4</a> ) ----> <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=4'>5</a> ctx.add_basemap(ax, <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=5'>6</a> crs=aoi.crs.to_string(), <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=6'>7</a> source=ctx.providers.CartoDB.Voyager <a href='vscode-notebook-cell:/c%3A/Users/jtrum/Desktop/maps/map_11_water-sources_children-elderly.ipynb#X42sZmlsZQ%3D%3D?line=7'>8</a> ) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\contextily\plotting.py:125, in add_basemap(ax, zoom, source, interpolation, attribution, attribution_size, reset_extent, crs, resampling, **extra_imshow_args) 123 # Convert extent from `crs` into WM for tile query 124 if crs is not None: --> 125 left, right, bottom, top = _reproj_bb( 126 left, right, bottom, top, crs, "epsg:3857" 127 ) 128 # Download image 129 image, extent = bounds2img( 130 left, bottom, right, top, zoom=zoom, source=source, ll=False 131 ) File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\contextily\plotting.py:215, in _reproj_bb(left, right, bottom, top, s_crs, t_crs) 214 def _reproj_bb(left, right, bottom, top, s_crs, t_crs): --> 215 n_l, n_b, n_r, n_t = transform_bounds(s_crs, t_crs, left, bottom, right, top) 216 return n_l, n_r, n_b, n_t File c:\Users\jtrum\miniconda3\envs\wash\lib\site-packages\rasterio\warp.py:147, in transform_bounds(src_crs, dst_crs, left, bottom, right, top, densify_pts) 113 def transform_bounds( 114 src_crs, 115 dst_crs, (...) 119 top, 120 densify_pts=21): 121 """Transform bounds from src_crs to dst_crs. 122 123 Optionally densifying the edges (to account for nonlinear transformations (...) 145 Outermost coordinates in target coordinate reference system. 146 """ --> 147 src_crs = CRS.from_user_input(src_crs) 148 dst_crs = CRS.from_user_input(dst_crs) 149 return _transform_bounds( 150 src_crs, 151 dst_crs, (...) 156 densify_pts, 157 ) File rasterio\\crs.pyx:777, in rasterio.crs.CRS.from_user_input() CRSError: The WKT could not be parsed. OGR Error code 6
In [ ]: