diff --git a/.nojekyll b/.nojekyll index faada17..367ff1d 100644 --- a/.nojekyll +++ b/.nojekyll @@ -1 +1 @@ -682ead62 \ No newline at end of file +00a753e4 \ No newline at end of file diff --git a/search.json b/search.json index c529096..3914528 100644 --- a/search.json +++ b/search.json @@ -130,7 +130,7 @@ "href": "tutorials/python/1-intro-python.html", "title": "Exploring the Legacy of Redlining", "section": "", - "text": "This executable notebook provides an opening example to illustrate a cloud-native workflow. Pedagogy research emphasizes the importance of “playing the whole game” before breaking down every pitch and hit. We intentionally focus on powerful high-level tools (STAC API, COGs, datacubes) to illustrate how a few chunks of code can perform a task that would be far slower and more verbose in a traditional file-based, download-first workflow. Note the close parallels between R and Python syntax. This arises because both languages wrap the same underlying tools (the STAC API and GDAL warper) and handle many of the nuisances of spatial data – from re-projections and resampling to mosaic tiles – without us noticing.\n\nfrom pystac_client import Client\nimport odc.stac\nimport pystac_client\nimport rioxarray\nimport os\n# suppress warning messages\nimport warnings\nwarnings.filterwarnings('ignore')\n\n\n\n\nThe first step in many workflows involves discovering individual spatial data files covering the space, time, and variables of interest. Here we use a STAC Catalog API to recover a list of candidate data. We dig deeper into how this works and what it returns in later recipes. This example searches for images in a lon-lat bounding box from a collection of Cloud-Optimized-GeoTIFF (COG) images taken by Sentinel2 satellite mission. This function will not download any imagery, it merely gives us a list of metadata about available images, including the access URLs.\n\nbox = [-122.51, 37.71, -122.36, 37.81]\nitems = (\n Client.\n open(\"https://earth-search.aws.element84.com/v1\").\n search(\n collections = ['sentinel-2-l2a'],\n bbox = box,\n datetime = \"2022-06-01/2022-08-01\",\n query={\"eo:cloud_cover\": {\"lt\": 20}}).\n item_collection()\n)\n\nWe pass this list of images to a high-level utilty (gdalcubes in R, odc.stac in python) that will do all of the heavy lifting. Using the URLs and metadata provided by STAC, these functions can extract only our data of interest (given by the bounding box) without downloading unnecessary regions or bands. While streaming the data, these functions will also reproject it into the desired coordinate reference system – (an often costly operation to perform in R) and can potentially resample or aggregate the data to a desired spatial resolution. (The R code will also resample from images in overlapping areas to replace pixels masked by clouds)\n\ndata = odc.stac.load(\n items,\n crs=\"EPSG:32610\",\n bands=[\"nir08\", \"red\"],\n bbox=box\n)\n# For some reason, requesting reprojection in EPSG:4326 gives only empty values\n\nWe can do arbitrary calculations on this data as well. Here we calculate NDVI, a widely used measure of greenness that can be used to determine vegetation cover. (Note that the R example uses lazy evaluation, and can thus perform these calculations while streaming)\n\nndvi = (((data.nir08 - data.red) / (data.red + data.nir08)).\n resample(time=\"MS\").\n median(\"time\", keep_attrs=True).\n compute()\n)\n\n# mask out bad pixels\nndvi = ndvi.where(ndvi <= 1)\n\n\nndvi\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n<xarray.DataArray (time: 1, y: 1118, x: 1329)> Size: 12MB\narray([[[nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n ...,\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan]]])\nCoordinates:\n * y (y) float64 9kB 4.185e+06 4.185e+06 ... 4.174e+06 4.174e+06\n * x (x) float64 11kB 5.431e+05 5.431e+05 ... 5.564e+05 5.564e+05\n spatial_ref int32 4B 32610\n * time (time) datetime64[ns] 8B 2022-06-01xarray.DataArraytime: 1y: 1118x: 1329nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nanarray([[[nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n ...,\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan]]])Coordinates: (4)y(y)float644.185e+06 4.185e+06 ... 4.174e+06units :metreresolution :-10.0crs :EPSG:32610array([4184925., 4184915., 4184905., ..., 4173775., 4173765., 4173755.])x(x)float645.431e+05 5.431e+05 ... 5.564e+05units :metreresolution :10.0crs :EPSG:32610array([543135., 543145., 543155., ..., 556395., 556405., 556415.])spatial_ref()int3232610spatial_ref :PROJCRS[\"WGS 84 / UTM zone 10N\",BASEGEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 10N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-123,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],BBOX[0,-126,84,-120]],ID[\"EPSG\",32610]]crs_wkt :PROJCRS[\"WGS 84 / UTM zone 10N\",BASEGEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 10N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-123,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],BBOX[0,-126,84,-120]],ID[\"EPSG\",32610]]semi_major_axis :6378137.0semi_minor_axis :6356752.314245179inverse_flattening :298.257223563reference_ellipsoid_name :WGS 84longitude_of_prime_meridian :0.0prime_meridian_name :Greenwichgeographic_crs_name :WGS 84horizontal_datum_name :World Geodetic System 1984 ensembleprojected_crs_name :WGS 84 / UTM zone 10Ngrid_mapping_name :transverse_mercatorlatitude_of_projection_origin :0.0longitude_of_central_meridian :-123.0false_easting :500000.0false_northing :0.0scale_factor_at_central_meridian :0.9996GeoTransform :543130 10 0 4184930 0 -10array(32610, dtype=int32)time(time)datetime64[ns]2022-06-01array(['2022-06-01T00:00:00.000000000'], dtype='datetime64[ns]')Indexes: (3)yPandasIndexPandasIndex(Index([4184925.0, 4184915.0, 4184905.0, 4184895.0, 4184885.0, 4184875.0,\n 4184865.0, 4184855.0, 4184845.0, 4184835.0,\n ...\n 4173845.0, 4173835.0, 4173825.0, 4173815.0, 4173805.0, 4173795.0,\n 4173785.0, 4173775.0, 4173765.0, 4173755.0],\n dtype='float64', name='y', length=1118))xPandasIndexPandasIndex(Index([543135.0, 543145.0, 543155.0, 543165.0, 543175.0, 543185.0, 543195.0,\n 543205.0, 543215.0, 543225.0,\n ...\n 556325.0, 556335.0, 556345.0, 556355.0, 556365.0, 556375.0, 556385.0,\n 556395.0, 556405.0, 556415.0],\n dtype='float64', name='x', length=1329))timePandasIndexPandasIndex(DatetimeIndex(['2022-06-01'], dtype='datetime64[ns]', name='time', freq='MS'))Attributes: (0)\n\n\nAnd we peek at the result 1. The long rectangle of Golden Gate Park is clearly visible in the North-West:\n\nndvi.plot()\n\n<matplotlib.collections.QuadMesh at 0x7c6ad6b50410>\n\n\n\n\n\nWe reproject this data into web mercator projection (longitude/latitude coordinates, EPSG:4326) and serialize as a Cloud Optimized Geotiff (COG) file using rioxarray, which will make this data useful to us in the future. While we could simply write the file to a path on the local disk (e.g. dest = \"sf-ndvi.tif\"), the code below exploits the virtual filesystem interface (VSI) to write the file directly to a cloud object store (which requires access to necessary environmental variables).\n\n# assumes configuration\nbucket = \"public-data\"\npath = \"espm-157/sf_ndvi.tif\"\ndest = f\"/vsis3/{bucket}/{path}\"\n\n( ndvi.\n rio.reproject(\"EPSG:4326\").\n rio.to_raster(dest, driver=\"COG\") \n)\n\nBecause we have written to a publicly readable bucket, our COG is now accessible without authentication using a normal https address. This will be convenient later in constructing interactive maps.\n\nendpoint = os.environ[\"AWS_S3_ENDPOINT\"]\npublic_url = f\"https://{endpoint}/{bucket}/{path}\"\n\n\n\n\n\nWe examine the present-day impact of historic “red-lining” of US cities during the Great Depression using data from the Mapping Inequality project. All though this racist practice was banned by federal law under the Fair Housing Act of 1968, the systemic scars of that practice are still so deeply etched on our landscape that the remain visible from space – “red-lined” areas (graded “D” under the racist HOLC scheme) show systematically lower greenness than predominately-white neighborhoods (Grade “A”). Trees provide many benefits, from mitigating urban heat to biodiversity, real-estate value, to health.\n\n\nWhile our satellite imagery data uses raster file formats such as COG, cartographic features such as polygons, lines or points are represented in vector formats. duckdb provides a scalable way to read and manipulate very large datasets from databases, tabular files (like csv or parquet) and geospatial vector objects (e.g. shapefiles, gpkg, geojson, geoparquet).\n\n\nimport ibis\nfrom ibis import _\ncon = ibis.duckdb.connect()\n\ncity = (\n con\n .read_geo(\"/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/mappinginequality.gpkg\")\n .filter(_.city == \"San Francisco\", _.residential)\n .execute()\n .set_crs(\"EPSG:4326\")\n )\n\n\n\n\n\n\n\nWe can use leafmap to visualize both raster and vector layers:\n\nimport leafmap.maplibregl as leafmap\n\nm = leafmap.Map() # set style=\"liberty\" for different basemap\nm.add_cog_layer(public_url, name = \"NDVI\", opacity = 0.7) # set palette=\"greens\" for color\nm.add_gdf(city, \"fill\", paint = {\"fill-color\": [\"get\", \"fill\"], \"fill-opacity\": 0.5})\nm.add_layer_control()\nm\n\n\n\n\n\n# use an iframe so the map can display on webpages\nfrom IPython.display import IFrame\nm.to_html(\"ndvi-redlines.html\", overwrite=True)\nIFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/tutorials/python/ndvi-redlines.html', width=700, height=400)\n\n\n \n \n\n\n\n\n\nIn addition to large scale raster data such as satellite imagery, the analysis of vector shapes such as polygons showing administrative regions is a central component of spatial analysis, and particularly important to spatial social sciences. The red-lined areas of the 1930s are one example of spatial vectors. One common operation is to summarise the values of all pixels falling within a given polygon, e.g. computing the average greenness (NDVI)\n\nfrom exactextract import exact_extract\nexact_extract(dest, \n city, \n [\"mean\"], \n include_geom = True,\n include_cols = [\"label\", \"grade\", \"city\", \"fill\"],\n output_options = {\"filename\": \"out.parquet\"},\n output=\"gdal\") # alternatively use \"pandas\" for (geo)pandas output\n\nAre historically redlined areas still less green?\n\n\nmean_ndvi = (con\n .read_parquet(\"out.parquet\")\n .rename(ndvi=\"mean\")\n .group_by(_.grade)\n .agg(mean_ndvi = _.ndvi.mean())\n .order_by(_.mean_ndvi.desc())\n .execute()\n)\n\n\nmean_ndvi\n\n\n\n\n\n\n\n\ngrade\nmean_ndvi\n\n\n\n\n0\nA\n0.338100\n\n\n1\nB\n0.247741\n\n\n2\nC\n0.236526\n\n\n3\nD\n0.232238\n\n\n\n\n\n\n\n\nimport altair as alt\n\n# grab the color scales from the data, for fun\ncolor_map = con.read_parquet(\"out.parquet\").select(\"grade\", \"fill\").distinct().order_by(_.grade).execute()\ncolors = alt.Scale(domain = color_map[\"grade\"], range = color_map[\"fill\"])\n\n# Create an Altair chart\nchart = alt.Chart(mean_ndvi).mark_bar().encode(\n x='grade',\n y='mean_ndvi:Q',\n color=alt.Color('grade', scale=colors)\n).properties(width=400)\n\n# display as SVG for compatibility with github ipynb rendering\nfrom IPython.display import SVG\nchart.save('temp.svg')\nSVG('temp.svg')" + "text": "This executable notebook provides an opening example to illustrate a cloud-native workflow. Pedagogy research emphasizes the importance of “playing the whole game” before breaking down every pitch and hit. We intentionally focus on powerful high-level tools (STAC API, COGs, datacubes) to illustrate how a few chunks of code can perform a task that would be far slower and more verbose in a traditional file-based, download-first workflow. Note the close parallels between R and Python syntax. This arises because both languages wrap the same underlying tools (the STAC API and GDAL warper) and handle many of the nuisances of spatial data – from re-projections and resampling to mosaic tiles – without us noticing.\n\nfrom pystac_client import Client\nimport odc.stac\nimport pystac_client\nimport rioxarray\nimport os\n# suppress warning messages\nimport warnings\nwarnings.filterwarnings('ignore')\n\n\n\n\nThe first step in many workflows involves discovering individual spatial data files covering the space, time, and variables of interest. Here we use a STAC Catalog API to recover a list of candidate data. We dig deeper into how this works and what it returns in later recipes. This example searches for images in a lon-lat bounding box from a collection of Cloud-Optimized-GeoTIFF (COG) images taken by Sentinel2 satellite mission. This function will not download any imagery, it merely gives us a list of metadata about available images, including the access URLs.\n\nbox = [-122.51, 37.71, -122.36, 37.81]\nitems = (\n Client.\n open(\"https://earth-search.aws.element84.com/v1\").\n search(\n collections = ['sentinel-2-l2a'],\n bbox = box,\n datetime = \"2022-06-01/2022-08-01\",\n query={\"eo:cloud_cover\": {\"lt\": 20}}).\n item_collection()\n)\n\nWe pass this list of images to a high-level utilty (gdalcubes in R, odc.stac in python) that will do all of the heavy lifting. Using the URLs and metadata provided by STAC, these functions can extract only our data of interest (given by the bounding box) without downloading unnecessary regions or bands. While streaming the data, these functions will also reproject it into the desired coordinate reference system – (an often costly operation to perform in R) and can potentially resample or aggregate the data to a desired spatial resolution. (The R code will also resample from images in overlapping areas to replace pixels masked by clouds)\n\ndata = odc.stac.load(\n items,\n crs=\"EPSG:32610\",\n bands=[\"nir08\", \"red\"],\n bbox=box\n)\n# For some reason, requesting reprojection in EPSG:4326 gives only empty values\n\nWe can do arbitrary calculations on this data as well. Here we calculate NDVI, a widely used measure of greenness that can be used to determine vegetation cover. (Note that the R example uses lazy evaluation, and can thus perform these calculations while streaming)\n\nndvi = (((data.nir08 - data.red) / (data.red + data.nir08)).\n resample(time=\"MS\").\n median(\"time\", keep_attrs=True).\n compute()\n)\n\n# mask out bad pixels\nndvi = ndvi.where(ndvi <= 1)\n\n\nndvi\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n<xarray.DataArray (time: 1, y: 1118, x: 1329)> Size: 12MB\narray([[[nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n ...,\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan]]])\nCoordinates:\n * y (y) float64 9kB 4.185e+06 4.185e+06 ... 4.174e+06 4.174e+06\n * x (x) float64 11kB 5.431e+05 5.431e+05 ... 5.564e+05 5.564e+05\n spatial_ref int32 4B 32610\n * time (time) datetime64[ns] 8B 2022-06-01xarray.DataArraytime: 1y: 1118x: 1329nan nan nan nan nan nan nan nan ... nan nan nan nan nan nan nan nanarray([[[nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n ...,\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan],\n [nan, nan, nan, ..., nan, nan, nan]]])Coordinates: (4)y(y)float644.185e+06 4.185e+06 ... 4.174e+06units :metreresolution :-10.0crs :EPSG:32610array([4184925., 4184915., 4184905., ..., 4173775., 4173765., 4173755.])x(x)float645.431e+05 5.431e+05 ... 5.564e+05units :metreresolution :10.0crs :EPSG:32610array([543135., 543145., 543155., ..., 556395., 556405., 556415.])spatial_ref()int3232610spatial_ref :PROJCRS[\"WGS 84 / UTM zone 10N\",BASEGEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 10N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-123,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],BBOX[0,-126,84,-120]],ID[\"EPSG\",32610]]crs_wkt :PROJCRS[\"WGS 84 / UTM zone 10N\",BASEGEOGCRS[\"WGS 84\",ENSEMBLE[\"World Geodetic System 1984 ensemble\",MEMBER[\"World Geodetic System 1984 (Transit)\"],MEMBER[\"World Geodetic System 1984 (G730)\"],MEMBER[\"World Geodetic System 1984 (G873)\"],MEMBER[\"World Geodetic System 1984 (G1150)\"],MEMBER[\"World Geodetic System 1984 (G1674)\"],MEMBER[\"World Geodetic System 1984 (G1762)\"],MEMBER[\"World Geodetic System 1984 (G2139)\"],MEMBER[\"World Geodetic System 1984 (G2296)\"],ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]],ENSEMBLEACCURACY[2.0]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433]],ID[\"EPSG\",4326]],CONVERSION[\"UTM zone 10N\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",-123,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",0,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]]],CS[Cartesian,2],AXIS[\"(E)\",east,ORDER[1],LENGTHUNIT[\"metre\",1]],AXIS[\"(N)\",north,ORDER[2],LENGTHUNIT[\"metre\",1]],USAGE[SCOPE[\"Navigation and medium accuracy spatial referencing.\"],AREA[\"Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK).\"],BBOX[0,-126,84,-120]],ID[\"EPSG\",32610]]semi_major_axis :6378137.0semi_minor_axis :6356752.314245179inverse_flattening :298.257223563reference_ellipsoid_name :WGS 84longitude_of_prime_meridian :0.0prime_meridian_name :Greenwichgeographic_crs_name :WGS 84horizontal_datum_name :World Geodetic System 1984 ensembleprojected_crs_name :WGS 84 / UTM zone 10Ngrid_mapping_name :transverse_mercatorlatitude_of_projection_origin :0.0longitude_of_central_meridian :-123.0false_easting :500000.0false_northing :0.0scale_factor_at_central_meridian :0.9996GeoTransform :543130 10 0 4184930 0 -10array(32610, dtype=int32)time(time)datetime64[ns]2022-06-01array(['2022-06-01T00:00:00.000000000'], dtype='datetime64[ns]')Indexes: (3)yPandasIndexPandasIndex(Index([4184925.0, 4184915.0, 4184905.0, 4184895.0, 4184885.0, 4184875.0,\n 4184865.0, 4184855.0, 4184845.0, 4184835.0,\n ...\n 4173845.0, 4173835.0, 4173825.0, 4173815.0, 4173805.0, 4173795.0,\n 4173785.0, 4173775.0, 4173765.0, 4173755.0],\n dtype='float64', name='y', length=1118))xPandasIndexPandasIndex(Index([543135.0, 543145.0, 543155.0, 543165.0, 543175.0, 543185.0, 543195.0,\n 543205.0, 543215.0, 543225.0,\n ...\n 556325.0, 556335.0, 556345.0, 556355.0, 556365.0, 556375.0, 556385.0,\n 556395.0, 556405.0, 556415.0],\n dtype='float64', name='x', length=1329))timePandasIndexPandasIndex(DatetimeIndex(['2022-06-01'], dtype='datetime64[ns]', name='time', freq='MS'))Attributes: (0)\n\n\nAnd we peek at the result 1. The long rectangle of Golden Gate Park is clearly visible in the North-West:\n\nndvi.plot()\n\n<matplotlib.collections.QuadMesh at 0x7c6ad6b50410>\n\n\n\n\n\nWe reproject this data into web mercator projection (longitude/latitude coordinates, EPSG:4326) and serialize as a Cloud Optimized Geotiff (COG) file using rioxarray, which will make this data useful to us in the future. While we could simply write the file to a path on the local disk (e.g. dest = \"sf-ndvi.tif\"), the code below exploits the virtual filesystem interface (VSI) to write the file directly to a cloud object store (which requires access to necessary environmental variables).\n\n# assumes configuration\nbucket = \"public-data\"\npath = \"espm-157/sf_ndvi.tif\"\ndest = f\"/vsis3/{bucket}/{path}\"\n\n( ndvi.\n rio.reproject(\"EPSG:4326\").\n rio.to_raster(dest, driver=\"COG\") \n)\n\nBecause we have written to a publicly readable bucket, our COG is now accessible without authentication using a normal https address. This will be convenient later in constructing interactive maps.\n\nendpoint = os.environ[\"AWS_S3_ENDPOINT\"]\npublic_url = f\"https://{endpoint}/{bucket}/{path}\"\n\n\n\n\n\nWe examine the present-day impact of historic “red-lining” of US cities during the Great Depression using data from the Mapping Inequality project. All though this racist practice was banned by federal law under the Fair Housing Act of 1968, the systemic scars of that practice are still so deeply etched on our landscape that the remain visible from space – “red-lined” areas (graded “D” under the racist HOLC scheme) show systematically lower greenness than predominately-white neighborhoods (Grade “A”). Trees provide many benefits, from mitigating urban heat to biodiversity, real-estate value, to health.\n\n\nWhile our satellite imagery data uses raster file formats such as COG, cartographic features such as polygons, lines or points are represented in vector formats. duckdb provides a scalable way to read and manipulate very large datasets from databases, tabular files (like csv or parquet) and geospatial vector objects (e.g. shapefiles, gpkg, geojson, geoparquet).\n\n\nimport ibis\nfrom ibis import _\ncon = ibis.duckdb.connect()\n\ncity = (\n con\n .read_geo(\"/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/mappinginequality.gpkg\")\n .filter(_.city == \"San Francisco\", _.residential)\n .execute()\n .set_crs(\"EPSG:4326\")\n )\n\n\n\n\n\n\n\nWe can use leafmap to visualize both raster and vector layers:\n\nimport leafmap.maplibregl as leafmap\n\nm = leafmap.Map() # set style=\"liberty\" for different basemap\nm.add_cog_layer(public_url, name = \"NDVI\", opacity = 0.7) # set palette=\"greens\" for color\nm.add_gdf(city, \"fill\", paint = {\"fill-color\": [\"get\", \"fill\"], \"fill-opacity\": 0.5})\nm.add_layer_control()\nm\n\n\n\n\n\n# use an iframe so the map can display on webpages\nfrom IPython.display import IFrame\nm.to_html(\"../../static/ndvi-redlines.html\", overwrite=True)\nIFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/ndvi-redlines.html', width=700, height=400)\n\n\n \n \n\n\n\n\n\nIn addition to large scale raster data such as satellite imagery, the analysis of vector shapes such as polygons showing administrative regions is a central component of spatial analysis, and particularly important to spatial social sciences. The red-lined areas of the 1930s are one example of spatial vectors. One common operation is to summarise the values of all pixels falling within a given polygon, e.g. computing the average greenness (NDVI)\n\nfrom exactextract import exact_extract\nexact_extract(dest, \n city, \n [\"mean\"], \n include_geom = True,\n include_cols = [\"label\", \"grade\", \"city\", \"fill\"],\n output_options = {\"filename\": \"out.parquet\"},\n output=\"gdal\") # alternatively use \"pandas\" for (geo)pandas output\n\nAre historically redlined areas still less green?\n\n\nmean_ndvi = (con\n .read_parquet(\"out.parquet\")\n .rename(ndvi=\"mean\")\n .group_by(_.grade)\n .agg(mean_ndvi = _.ndvi.mean())\n .order_by(_.mean_ndvi.desc())\n .execute()\n)\n\n\nmean_ndvi\n\n\n\n\n\n\n\n\ngrade\nmean_ndvi\n\n\n\n\n0\nA\n0.338100\n\n\n1\nB\n0.247741\n\n\n2\nC\n0.236526\n\n\n3\nD\n0.232238\n\n\n\n\n\n\n\n\nimport altair as alt\n\n# grab the color scales from the data, for fun\ncolor_map = con.read_parquet(\"out.parquet\").select(\"grade\", \"fill\").distinct().order_by(_.grade).execute()\ncolors = alt.Scale(domain = color_map[\"grade\"], range = color_map[\"fill\"])\n\n# Create an Altair chart\nchart = alt.Chart(mean_ndvi).mark_bar().encode(\n x='grade',\n y='mean_ndvi:Q',\n color=alt.Color('grade', scale=colors)\n).properties(width=400)\n\n# display as SVG for compatibility with github ipynb rendering\nfrom IPython.display import SVG\nchart.save('temp.svg')\nSVG('temp.svg')" }, { "objectID": "tutorials/python/1-intro-python.html#data-discovery", @@ -151,7 +151,7 @@ "href": "tutorials/python/1-intro-python.html#interactive-maps", "title": "Exploring the Legacy of Redlining", "section": "", - "text": "We can use leafmap to visualize both raster and vector layers:\n\nimport leafmap.maplibregl as leafmap\n\nm = leafmap.Map() # set style=\"liberty\" for different basemap\nm.add_cog_layer(public_url, name = \"NDVI\", opacity = 0.7) # set palette=\"greens\" for color\nm.add_gdf(city, \"fill\", paint = {\"fill-color\": [\"get\", \"fill\"], \"fill-opacity\": 0.5})\nm.add_layer_control()\nm\n\n\n\n\n\n# use an iframe so the map can display on webpages\nfrom IPython.display import IFrame\nm.to_html(\"ndvi-redlines.html\", overwrite=True)\nIFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/tutorials/python/ndvi-redlines.html', width=700, height=400)" + "text": "We can use leafmap to visualize both raster and vector layers:\n\nimport leafmap.maplibregl as leafmap\n\nm = leafmap.Map() # set style=\"liberty\" for different basemap\nm.add_cog_layer(public_url, name = \"NDVI\", opacity = 0.7) # set palette=\"greens\" for color\nm.add_gdf(city, \"fill\", paint = {\"fill-color\": [\"get\", \"fill\"], \"fill-opacity\": 0.5})\nm.add_layer_control()\nm\n\n\n\n\n\n# use an iframe so the map can display on webpages\nfrom IPython.display import IFrame\nm.to_html(\"../../static/ndvi-redlines.html\", overwrite=True)\nIFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/ndvi-redlines.html', width=700, height=400)" }, { "objectID": "tutorials/python/1-intro-python.html#zonal-statistics", diff --git a/sitemap.xml b/sitemap.xml index 3393a6d..72ebad5 100644 --- a/sitemap.xml +++ b/sitemap.xml @@ -2,26 +2,26 @@ https://github.com/boettiger-lab/nasa-topst-env-justice/tutorials/python/2-NASA-EarthData.html - 2024-12-29T04:37:12.820Z + 2024-12-29T04:51:25.609Z https://github.com/boettiger-lab/nasa-topst-env-justice/tutorials/computing-environment.html - 2024-12-29T04:37:11.880Z + 2024-12-29T04:51:24.674Z https://github.com/boettiger-lab/nasa-topst-env-justice/tutorials/R/1-intro-R.html - 2024-12-29T04:37:10.946Z + 2024-12-29T04:51:23.773Z https://github.com/boettiger-lab/nasa-topst-env-justice/index.html - 2024-12-29T04:37:10.080Z + 2024-12-29T04:51:22.918Z https://github.com/boettiger-lab/nasa-topst-env-justice/tutorials/R/2-earthdata.html - 2024-12-29T04:37:11.526Z + 2024-12-29T04:51:24.330Z https://github.com/boettiger-lab/nasa-topst-env-justice/tutorials/python/1-intro-python.html - 2024-12-29T04:37:12.337Z + 2024-12-29T04:51:25.129Z diff --git a/static/ndvi-redlines.html b/static/ndvi-redlines.html new file mode 100644 index 0000000..290a9a2 --- /dev/null +++ b/static/ndvi-redlines.html @@ -0,0 +1,38 @@ + + + +My Awesome Map + + + + + + +
+ + + \ No newline at end of file diff --git a/tutorials/python/1-intro-python.html b/tutorials/python/1-intro-python.html index 785189e..e065ee4 100644 --- a/tutorials/python/1-intro-python.html +++ b/tutorials/python/1-intro-python.html @@ -863,14 +863,14 @@

Interactive maps

-
+
# use an iframe so the map can display on webpages
 from IPython.display import IFrame
-m.to_html("ndvi-redlines.html", overwrite=True)
-IFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/tutorials/python/ndvi-redlines.html', width=700, height=400)
-
+m.to_html("../../static/ndvi-redlines.html", overwrite=True) +IFrame(src='https://boettiger-lab.github.io/nasa-topst-env-justice/ndvi-redlines.html', width=700, height=400)
+
- +