wip; tinkering hyperparameters for large geoms

preprocess/CBN-data.ipynb

@@ -63,11 +63,11 @@
     "folder = 'Counties'\n",
     "name = 'CA_counties'\n",
     "\n",
-    "# unzip(folder = folder, file = '30x30_Counties.zip')\n",
-    "# process_vector(folder = folder, file = f\"{name}.shp\")\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.parquet\")\n",
+    "unzip(s3, folder = folder, file = '30x30_Counties.zip')\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= \"COUNTY_NAM\")"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)"
    ]
   },
   {
@@ -91,9 +91,9 @@
     "\n",
     "folder = 'Climate_zones'\n",
     "name = 'climate_zones_10'\n",
-    "# download(folder = folder, file = 'clusters_10.tif')\n",
-    "# process_raster(s3, folder = folder, file = 'clusters_10.tif', file_name = f\"{name}.tif\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols= \"id\")"
+    "download(folder = folder, file = 'clusters_10.tif')\n",
+    "cols = process_raster(s3, folder = folder, file = 'clusters_10.tif', file_name = f\"{name}.tif\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols)"
    ]
   },
   {
@@ -118,10 +118,10 @@
     "folder = 'Ecoregion'\n",
     "name = 'ACE_ecoregions'\n",
     "\n",
-    "unzip(folder = folder, file = '30x30_Ecoregions.zip')\n",
-    "process_vector(folder = folder, file = f\"{name}.shp\")\n",
+    "unzip(s3, folder = folder, file = '30x30_Ecoregions.zip')\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\")\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= 'CA_Ecoregi')"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)"
    ]
   },
   {
@@ -148,7 +148,7 @@
    "outputs": [],
    "source": [
     "# download(folder = 'Habitat', file = 'CWHR13_2022.tif')\n",
-    "# process_raster(s3, folder = 'Habitat', file = 'CWHR13_2022.tif')"
+    "# cols = process_raster(s3, folder = 'Habitat', file = 'CWHR13_2022.tif')"
    ]
   },
   {
@@ -165,20 +165,20 @@
     "folder = 'Habitat'\n",
     "name = 'fveg22_1'\n",
     "\n",
-    "# unzip(folder = folder, file = 'fveg221gdb.zip')\n",
+    "unzip(s3, folder = folder, file = 'fveg221gdb.zip')\n",
     "\n",
-    "# command = [\n",
-    "#     \"gdalwarp\",\n",
-    "#     \"-of\", \"GTiff\",\n",
-    "#     'fveg22_1.gdb',\n",
-    "#     'fveg22_1.tif' \n",
-    "#     ]\n",
+    "command = [\n",
+    "    \"gdalwarp\",\n",
+    "    \"-of\", \"GTiff\",\n",
+    "    'fveg22_1.gdb',\n",
+    "    'fveg22_1.tif' \n",
+    "    ]\n",
     "\n",
-    "# subprocess.run(command, check=True)\n",
-    "process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "subprocess.run(command, check=True)\n",
+    "cols = process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
     "upload(folder = folder, file = f'{name}_processed.tif.aux.xml')\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols= \"id\")"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols)"
    ]
   },
   {
@@ -211,11 +211,11 @@
     "download(folder = folder, file = 'Terrestrial_Biodiversity_Summary_-_ACE_[ds2739].geojson',\n",
     "         file_name = f\"{name}.geojson\")\n",
     "\n",
-    "process_vector(folder = folder, file = f\"{name}.geojson\")\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.geojson\")\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.geojson\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.geojson\")\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= \"OBJECTID\")\n",
-    "# gdf = gpd.read_parquet(f\"{name}.parquet\")\n"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n",
+    "gdf = gpd.read_parquet(f\"{name}.parquet\")\n"
    ]
   },
   {
@@ -244,8 +244,8 @@
     "            'County', 'Shape__Area', 'Shape__Length', 'geometry']\n",
     "    cols.append(col) #select only the cols we want + the new col. \n",
     "    rank_df = gdf[gdf[col]==5][cols]# filter ranks = 5\n",
-    "    process_vector(folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = rank_df)\n",
-    "    convert_pmtiles(folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n"
+    "    cols = process_vector(s3, folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = rank_df)\n",
+    "    convert_pmtiles(con, s3, folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n"
    ]
   },
   {
@@ -283,14 +283,32 @@
     "        percentile = 0.95\n",
     "    threshold = gdf[col].quantile(percentile)\n",
     "    ace = gdf[gdf[col]>=threshold][cols]\n",
-    "    process_vector(folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = ace)\n",
-    "    convert_pmtiles(folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n",
+    "    cols = process_vector(s3, folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = ace)\n",
+    "    convert_pmtiles(con, s3, folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n",
     "\n",
     "\n",
     "# calculate 80% percentile, filter to those >= threshold. \n",
     "# subset to calculate acres within each network, % of feature conserved and % of network "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50f9c3bc-8e7e-4bb9-b1c9-9718cf8454a8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "con = ibis.duckdb.connect(extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "folder = 'Climate_risks/Historical_fire_perimeters'\n",
+    "name = 'calfire_2023'\n",
+    "\n",
+    "url = f\"s3://public-ca30x30/CBN-data/{folder}/hex/{name}.parquet\"\n",
+    "\n",
+    "con.read_parquet(url).head(10).execute()\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "6991222f-7d24-4f10-9ee0-db20513405d6",
@@ -321,9 +339,9 @@
     "folder = 'Biodiversity_unique/Plant_richness'\n",
     "name = 'species_D'\n",
     "\n",
-    "download(folder = folder, file = f\"{name}.tif\")\n",
-    "filter_raster(folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols= \"id\")"
+    "# download(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "cols = filter_raster(s3, folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols)"
    ]
   },
   {
@@ -348,9 +366,9 @@
     "folder = 'Biodiversity_unique/Rarityweighted_endemic_plant_richness'\n",
     "name = 'endemicspecies_E'\n",
     "\n",
-    "download(folder = folder, file = f\"{name}.tif\")\n",
-    "filter_raster(folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols= \"id\")"
+    "download(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "cols = filter_raster(s3, folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols)"
    ]
   },
   {
@@ -392,8 +410,8 @@
     "folder = 'Connectivity_resilience/Resilient_connected_network_allcategories'\n",
     "name = 'rcn_wIntactBioCat_caOnly_2020-10-27'\n",
     "\n",
-    "process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols= \"id\")"
+    "cols = process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols)"
    ]
   },
   {
@@ -468,20 +486,21 @@
    "outputs": [],
    "source": [
     "%%time \n",
+    "con = ibis.duckdb.connect('wetlands',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
     "\n",
     "folder = 'Freshwater_resources/Wetlands'\n",
     "name = 'CA_wetlands'\n",
     "\n",
     "# only pick a subset \n",
-    "unzip(folder = folder, file = 'CA_geodatabase_wetlands.zip')\n",
+    "unzip(s3, folder = folder, file = 'CA_geodatabase_wetlands.zip')\n",
     "gdf = gpd.read_file('CA_geodatabase_wetlands.gdb')\n",
     "wetlands = ['Freshwater Emergent Wetland', 'Freshwater Forested/Shrub Wetland', 'Estuarine and Marine Wetland']\n",
     "gdf = gdf[gdf['WETLAND_TYPE'].isin(wetlands)]\n",
     "\n",
-    "process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
-    "# convert_pmtiles(folder =folder, file = f\"{name}.parquet\")\n",
-    "geom_to_h3(con, folder = folder, file = f\"{name}.parquet\", cols= ['ATTRIBUTE','WETLAND_TYPE','NWI_ID'])\n",
-    "\n"
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
+    "convert_pmtiles(con, s3, folder =folder, file = f\"{name}.parquet\")\n",
+    "geom_to_h3(con, folder = folder, file = f\"{name}.parquet\", cols = cols)\n"
    ]
   },
   {
@@ -580,34 +599,33 @@
    "outputs": [],
    "source": [
     "%%time \n",
+    "con = ibis.duckdb.connect('farm',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
     "\n",
     "folder = 'NBS_agriculture/Farmland'\n",
-    "unzip(folder = folder, file = 'Important_Farmland_2018.zip')\n",
+    "unzip(s3, folder = folder, file = 'Important_Farmland_2018.zip')\n",
     "\n",
     "folder = 'NBS_agriculture/Farmland_all'\n",
     "name = 'Important_Farmland_2018'\n",
-    "process_vector(folder = folder, file = f\"{name}.gdb\")\n",
-    "# convert_pmtiles(folder = folder, file =f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= ['county_nam','polygon_ty'])\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.gdb\",crs = \"epsg:4326\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n",
     "\n",
     "# only pick a subset \n",
     "folder = 'NBS_agriculture/Farmland_all/Farmland'\n",
     "name = 'Farmland_2018'\n",
-    "# gdf = gpd.read_file('Important_Farmland_2018.gdb')\n",
-    "# farmland_type = ['P','S','L','U'] # prime, statewide importance, local importance, unique\n",
-    "# gdf_farmland = gdf[gdf['polygon_ty'].isin(farmland_type)]\n",
-    "# process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf_farmland)\n",
-    "# convert_pmtiles(folder = folder, file =f\"{name}.parquet\")\n",
-    "\n",
-    "\n",
+    "gdf = gpd.read_file('Important_Farmland_2018.gdb')\n",
+    "farmland_type = ['P','S','L','U'] # prime, statewide importance, local importance, unique\n",
+    "gdf_farmland = gdf[gdf['polygon_ty'].isin(farmland_type)]\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf_farmland)\n",
+    "convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n",
     "\n",
     "# grazing lands \n",
     "folder = 'NBS_agriculture/Farmland_all/Lands_suitable_grazing'\n",
     "name = 'Grazing_land_2018'\n",
-    "\n",
-    "# gdf_grazing = gdf[gdf['polygon_ty'] == 'G']\n",
-    "# process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf_grazing)\n",
-    "# convert_pmtiles(folder = folder, file =f\"{name}.parquet\")\n"
+    "gdf_grazing = gdf[gdf['polygon_ty'] == 'G']\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf_grazing)\n",
+    "convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n"
    ]
   },
   {
@@ -644,14 +662,6 @@
     "Only YEAR >= 2014. "
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "425f9149-d8ac-437a-9572-301bd1b1bec8",
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -663,15 +673,15 @@
     "folder = 'Climate_risks/Historical_fire_perimeters'\n",
     "name = 'calfire_2023'\n",
     "\n",
-    "unzip(folder = folder, file = 'fire23-1gdb.zip')\n",
+    "unzip(s3, folder = folder, file = 'fire23-1gdb.zip')\n",
     "gdf = gpd.read_file('fire23_1.gdb')\n",
-    "gdf = gdf[~gdf['YEAR_'].isna()]\n",
-    "gdf['YEAR_'] = gdf['YEAR_'].astype('int64')\n",
+    "# gdf = gdf[~gdf['YEAR_'].isna()]\n",
+    "# gdf['YEAR_'] = gdf['YEAR_'].astype('int64')\n",
     "# gdf = gdf[gdf['YEAR_']>=2014]\n",
-    "process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.parquet\")\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= ['INC_NUM','FIRE_NAME','YEAR_'])"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)"
    ]
   },
   {
@@ -720,7 +730,7 @@
     "Do seperately for both climate models - CNRM and MIROC.\n",
     "'''\n",
     "\n",
-    "unzip(folder = 'Climate_risks/Mid-century_habitat_climate_exposure', file = 'Midcentury_habitat_climate_exposure.zip')\n",
+    "unzip(s3, folder = 'Climate_risks/Mid-century_habitat_climate_exposure', file = 'Midcentury_habitat_climate_exposure.zip')\n",
     "\n",
     "# still need to do "
    ]
@@ -752,10 +762,10 @@
     "folder = 'Progress_data_new_protection/Newly_counted_lands'\n",
     "name = 'newly_counted_lands_2024'\n",
     "\n",
-    "unzip(folder = folder, file = f\"{name}.shp.zip\")\n",
-    "process_vector(folder = folder, file = f\"{name}.shp\")\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= ['ORIG_FID'])\n"
+    "unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n"
    ]
   },
   {
@@ -777,10 +787,10 @@
     "folder = 'Progress_data_new_protection/DAC'\n",
     "name = 'DAC_2022'\n",
     "\n",
-    "unzip(folder = folder, file = 'sb535dacgdbf2022gdb.zip')\n",
-    "process_vector(folder = folder, file = 'SB535DACgdb_F_2022.gdb', file_name = f\"{name}.parquet\")\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= ['GEOID'])\n"
+    "unzip(s3, folder = folder, file = 'sb535dacgdbf2022gdb.zip')\n",
+    "cols = process_vector(s3, folder = folder, file = 'SB535DACgdb_F_2022.gdb', file_name = f\"{name}.parquet\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n"
    ]
   },
   {
@@ -804,34 +814,17 @@
     "\n",
     "folder = 'Progress_data_new_protection/Priority_populations'\n",
     "name = 'CalEnviroScreen4'\n",
-    "unzip(folder = folder, file = 'Priority Populations 4.0 Geodatabase.zip')\n",
+    "unzip(s3, folder = folder, file = 'Priority Populations 4.0 Geodatabase.zip')\n",
     "\n",
     "gdf = (con.read_geo('Priority Populations 4.0 Combined Layer.gdb')\n",
     "       .mutate(id=ibis.row_number().over()) #making a unique id \n",
     "      ).execute().set_crs('EPSG:3857')\n",
     "\n",
-    "process_vector(folder = folder, file = 'Priority Populations 4.0 Combined Layer.gdb',\n",
+    "cols = process_vector(folder = folder, file = 'Priority Populations 4.0 Combined Layer.gdb',\n",
     "               file_name = f\"{name}.parquet\", gdf = gdf)\n",
     "\n",
-    "# convert_pmtiles(folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= [\"id\"])\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "e64129da-f369-425f-afcc-bc595a89fb7d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "file = f\"{name}.parquet\"\n",
-    "folder = 'Progress_data_new_protection/Priority_populations'\n",
-    "name = 'CalEnviroScreen4'\n",
-    "bucket, path = info(folder, file)\n",
-    "# path, file = os.path.split(path)\n",
-    "# name, ext = os.path.splitext(file)\n",
-    "# join_chunked(bucket, path, file)\n",
-    "con.read_parquet(f\"s3://{bucket}/{folder}/hex/{file}_part_000.parquet\").head(10).execute()"
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n"
    ]
   },
   {
@@ -852,12 +845,12 @@
     "folder = 'Progress_data_new_protection/Low_income_communities'\n",
     "name = 'low_income_CalEnviroScreen4'\n",
     "\n",
-    "unzip(folder = folder, file = 'Priority Populations 4.0 Geodatabase.zip')\n",
+    "unzip(s3, folder = folder, file = 'Priority Populations 4.0 Geodatabase.zip')\n",
     "\n",
     "gdf = gpd.read_file('Priority Populations 4.0 Combined Layer.gdb')\n",
     "gdf = gdf[gdf['Designatio'] =='Low-income community']\n",
-    "process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
-    "convert_pmtiles(folder = folder, file = f\"{name}.parquet\")"
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")"
    ]
   },
   {
@@ -878,10 +871,10 @@
     "folder = 'Progress_data_new_protection/Land_Status_Zone_Ecoregion_Counties'\n",
     "name = 'all_regions_reGAP_county_eco'\n",
     "\n",
-    "unzip(folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp.zip')\n",
-    "process_vector(folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp',\n",
+    "unzip(s3, folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp.zip')\n",
+    "cols = process_vector(s3, folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp',\n",
     "               file_name = f\"{name}.parquet\")\n",
-    "convert_pmtiles(folder = folder, file = f\"{name}.parquet\")"
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")"
    ]
   },
   {
@@ -892,43 +885,6 @@
     "# CA Nature data"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "ecc0f168-badd-4e4d-b97b-ee7891afaa4e",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# def convert_h3_2(con, folder, file, cols, zoom = \"8\"):\n",
-    "#     cols = \", \".join(cols) if isinstance(cols,list) else cols #unpack columns  \n",
-    "#     bucket = 'public-ca30x30'\n",
-    "#     name = 'ca-30x30-base'\n",
-    "#     file = 'ca-30x30-base.parquet'\n",
-    "#     folder = \"Preprocessing\"\n",
-    "#     name= name.replace('-','')\n",
-    "\n",
-    "\n",
-    "#     # reproject \n",
-    "#     # (con.read_parquet(f\"s3://{bucket}/{file}\")\n",
-    "#     #  .mutate(geom = _.geom.convert('epsg:3310','epsg:4326'))\n",
-    "#     # ).to_parquet(f\"s3://{bucket}/hex/{file}\")\n",
-    "\n",
-    "#     con.read_parquet(f\"s3://{bucket}/{folder}/{file}\", table_name = name)\n",
-    "\n",
-    "#     con.sql(f'''\n",
-    "#       WITH t2 AS (\n",
-    "#         WITH t1 AS (\n",
-    "#           SELECT {cols},  ST_Dump(geom) AS geom \n",
-    "#           FROM {name}\n",
-    "#         ) \n",
-    "#         SELECT {cols},\n",
-    "#               h3_polygon_wkt_to_cells_string(UNNEST(geom).geom, {zoom}) AS h{zoom}\n",
-    "#         FROM t1\n",
-    "#       )\n",
-    "#       SELECT *, UNNEST(h{zoom}) AS h{zoom} FROM t2\n",
-    "#       ''').to_parquet(f\"s3://{bucket}/{folder}/hex/{file}\")"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -946,46 +902,44 @@
     "# download(folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
     "# gdf = gpd.read_parquet(f\"{name}.parquet\")\n",
-    "process_vector(folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= [\"id\"])\n",
-    "\n",
-    "# convert_h3_2(con, folder = folder, file = f\"{name}.parquet\", cols= [\"id\"])\n",
-    "\n"
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols)\n"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "908786ec-2a86-4fb0-a47a-9de364254806",
+   "cell_type": "markdown",
+   "id": "af486c71-3b84-4685-9794-fbacbf5f81c7",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# url = f\"s3://public-ca30x30/{folder}/hex/{name}.parquet\"\n",
-    "# # url = f\"s3://public-ca30x30/{folder}/{name}.parquet\"\n",
-    "\n",
-    "# con.read_parquet(url).head(10).execute()\n"
+    "# CPAD"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "65d98aa3-041f-42d6-8448-6fa8a05f5850",
+   "id": "cf6c896f-65f3-403a-abd9-f7dec2f4f112",
    "metadata": {},
    "outputs": [],
    "source": [
-    "# file = 'ca-30x30-base.parquet'\n",
-    "# folder = \"Preprocessing\"\n",
-    "# bucket = 'public-ca30x30'\n",
-    "# con.read_parquet(f\"s3://{bucket}/{folder}/hex/{file}\").head(10).execute()"
+    "con = ibis.duckdb.connect('cpad',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "folder = 'cpad'\n",
+    "name = 'cced_2024b_release'\n",
+    "\n",
+    "# unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "# cols = process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:3310\")\n",
+    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:4326\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= cols)\n",
+    "\n",
+    "name = 'cpad_2024b_release'\n",
+    "# unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "# cols = process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:3310\")\n",
+    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:4326\")\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= cols)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "fa960a99-3c79-4e67-becf-a1cb397aa5fb",
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {

preprocess/h3_utils.py

@@ -1,5 +1,16 @@
+from utils import *
 
-def geom_to_h3(con, name, cols, zoom):
+# === CONFIG ===
+default_zoom = "8"
+default_limit = 10_000
+default_geom_len_thresh = 5_000  # H3 cells per geometry
+chunk_limit = default_limit
+large_geom_thresh = default_geom_len_thresh
+est_total_h3_thresh = 150_000 
+large_geom_batch_limit = 100
+
+
+def compute_h3(con, name, cols, zoom):
     """
     Computes hexes    
     """
@@ -38,38 +49,36 @@ def check_size(con, name, zoom, sample_size=100):
 
     return est_total_h3, max_len
 
+# def chunk_large_geom(con, s3, bucket, path, name, zoom=default_zoom, geom_len_threshold=large_geom_thresh):
+# def chunk_large_geom(con, s3, bucket, path, name, zoom="8", geom_len_threshold=10_000):
 
-def write_large_geoms(con, s3, bucket, path, name, zoom="8", geom_len_threshold=10_000):
+def chunk_large_geom(con, s3, bucket, path, name, zoom=default_zoom,
+                      geom_len_threshold=large_geom_thresh,
+                      batch_limit=large_geom_batch_limit):
     """
     Individually processing large geoms (different from processing "chunks")
     """
     offset = 0
     i = 0
-    limit=3000
     while True:
-        large_key = f"{path}/hex/{name}_large_{i:03d}.parquet"
-        print(f"🟠 Checking large geometry batch {i} → {large_key}")
+        relative_key = f"{path}/hex/{name}_large_{i:03d}.parquet"
+        print(f"🟠 Checking large geometry batch {i} → {relative_key}")
 
-        #  check if file already exists in minio
-        try:
-            s3.stat_object(bucket, large_key)
-            print(f"⏩ Skipping existing large batch: {large_key}")
-            offset += limit
+        if exists_on_s3(s3, folder="", file=relative_key):  # we pass relative_key as `file`
+            print(f"⏩ Skipping existing large batch: {relative_key}")
+            offset += batch_limit
             i += 1
             continue
-        except S3Error as err:
-            if err.code != "NoSuchKey":
-                raise  
 
-        print(f"📝 Writing large geometry batch {i} → {large_key}")
+        print(f"📝 Writing large geometry batch {i} → {relative_key}")
         q = con.sql(f'''
             SELECT *, UNNEST(h{zoom}) AS h{zoom}
             FROM t2
             WHERE len(h{zoom}) > {geom_len_threshold}
-            LIMIT {limit} OFFSET {offset}
+            LIMIT {batch_limit} OFFSET {offset}
         ''')
 
-        q.to_parquet(f"s3://{bucket}/{large_key}")
+        q.to_parquet(f"s3://{bucket}/{relative_key}")
 
         if q.count().execute() == 0:
             break
@@ -79,7 +88,6 @@ def write_large_geoms(con, s3, bucket, path, name, zoom="8", geom_len_threshold=
 
     return i
 
-    
 def join_large_geoms(con, s3, bucket, path, name):
     """
     If we had to process large geoms individually, join those datasets after conversion.
@@ -87,14 +95,10 @@ def join_large_geoms(con, s3, bucket, path, name):
     # check if any large files exist before trying to join
     test_key = f"{path}/hex/{name}_large_000.parquet"
 
-    try:
-        s3.stat_object(bucket, test_key)
-    except S3Error as err:
-        if err.code == "NoSuchKey":
-            print("✅ No large geometry chunks to join.")
-            return
-        else:
-            raise
+    if not exists_on_s3(s3, folder="", file=test_key):
+        print("✅ No large geometry chunks to join.")
+        return
+
     # join if it exists 
     con.raw_sql(f'''
         COPY (
@@ -103,27 +107,24 @@ def join_large_geoms(con, s3, bucket, path, name):
         TO 's3://{bucket}/{path}/hex/{name}_large.parquet'
         (FORMAT PARQUET)
     ''')
+    
 
-
-def chunk_data(con, s3, bucket, path, name, zoom="8", limit=100_000, geom_len_threshold=10_000):
+# def chunk_geom(con, s3, bucket, path, name, zoom="8", limit=50_000, geom_len_threshold=10_000):
+def chunk_geom(con, s3, bucket, path, name, zoom=default_zoom, limit=chunk_limit, geom_len_threshold=large_geom_thresh):
     """
     Processing large files in chunks. 
     """
     offset = 0
     i = 0
-
+    
     while True:
         chunk_path = f"{path}/hex/{name}_chunk{i:03d}.parquet"
         
-        try:
-            s3.stat_object(bucket, chunk_path)
+        if exists_on_s3(s3, folder="", file=chunk_path):  # relative path passed as file
             print(f"⏩ Skipping existing chunk: {chunk_path}")
             offset += limit
             i += 1
             continue
-        except S3Error as err:
-            if err.code != "NoSuchKey":
-                raise
 
         print(f"📝 Writing chunk {i} → {chunk_path}")
         q = con.sql(f'''
@@ -139,13 +140,13 @@ def chunk_data(con, s3, bucket, path, name, zoom="8", limit=100_000, geom_len_th
         i += 1
 
     # process large geometries using same threshold and limit
-    write_large_geoms(con, s3, bucket, path, name, zoom, geom_len_threshold=geom_len_threshold)
+    chunk_large_geom(con, s3, bucket, path, name, zoom, geom_len_threshold=geom_len_threshold)
     join_large_geoms(con, s3, bucket, path, name)
     return i
 
 
 
-def join_chunked(bucket, path, name):
+def join_chunked(con, bucket, path, name):
     """
     If we had to chunk the data, join those datasets after conversion.
     """
@@ -158,7 +159,8 @@ def join_chunked(bucket, path, name):
         ''')
 
 # def convert_h3(con, folder, file, cols, zoom="8", limit=100_000, geom_len_threshold=10_000):
-def convert_h3(con, s3, folder, file, cols, zoom="8", limit=100_000, geom_len_threshold=5_000):
+# def convert_h3(con, s3, folder, file, cols, zoom="8", limit=50_000, geom_len_threshold=5_000):
+def convert_h3(con, s3, folder, file, cols, zoom=default_zoom, limit=chunk_limit, geom_len_threshold=large_geom_thresh):
     """
     Driver function to convert geometries to h3
     """
@@ -175,14 +177,14 @@ def convert_h3(con, s3, folder, file, cols, zoom="8", limit=100_000, geom_len_th
     est_total, max_per_geom = check_size(con, name, zoom)
     # if est_total > 500_000 or max_per_geom > geom_len_threshold:
 
-    if est_total > 1_000_000 or max_per_geom > geom_len_threshold:
+    if est_total > est_total_h3_thresh or max_per_geom > geom_len_threshold:
         print("Chunking due to estimated size")
-        geom_to_h3(con, name, cols, zoom)
-        chunk_data(con, s3, bucket, path, name, zoom, limit, geom_len_threshold)
+        compute_h3(con, name, cols, zoom)
+        chunk_geom(con, s3, bucket, path, name, zoom, limit, geom_len_threshold)
         join_chunked(con, bucket, path, name)
     else:
         print("Writing single output")
-        geom_to_h3(con, name, cols, zoom)
+        compute_h3(con, name, cols, zoom)
         con.sql(f'''
             SELECT *, UNNEST(h{zoom}) AS h{zoom}
             FROM t2

preprocess/utils.py

@@ -1,4 +1,5 @@
 from minio.error import S3Error
+from cng.utils import *
 
 import zipfile
 import os
@@ -40,7 +41,7 @@ def upload(s3, folder, file):
     s3.fput_object(bucket, path ,file) 
     return
 
-def unzip(folder, file):
+def unzip(s3, folder, file):
     """
     Unzipping zip files 
     """
@@ -49,8 +50,8 @@ def unzip(folder, file):
         zip_ref.extractall()
     return 
 
-# def process_vector(folder, file, file_name = None, gdf = None, crs="EPSG:3310"):
-def process_vector(folder, file, file_name = None, gdf = None, crs="EPSG:4326"):
+# def process_vector(s3, folder, file, file_name = None, gdf = None, crs="EPSG:3310"):
+def process_vector(s3, folder, file, file_name = None, gdf = None, crs="EPSG:4326"):
     """
     Driver function to process vectors 
     """
@@ -67,20 +68,56 @@ def process_vector(folder, file, file_name = None, gdf = None, crs="EPSG:4326"):
     parquet_file = f"{name}{'.parquet'}"
     gdf.to_parquet(parquet_file)
     upload(s3, folder, parquet_file)
-    return 
-
 
-# def upload_parquet(folder, file, gdf):
-#     """
-#     Uploading parquets 
-#     """
-#     name, ext = os.path.splitext(file)
-#     parquet_file = f"{name}{'.parquet'}"
-#     gdf.to_parquet(parquet_file)
-#     upload(folder, parquet_file)
-#     return  
+    return gdf.drop('geom',axis = 1).columns.to_list()
 
+def process_raster(s3, folder, file, file_name = None):
+    """
+    Driver function to process rasters 
+    """
+    if file_name:
+        file = file_name
+    # output_file = reproject_raster(file)
+    # upload(s3, folder, output_file)
+    # output_cog_file = make_cog(output_file)
+    # upload(s3, folder, output_cog_file)
+    # output_vector, cols  = make_vector(output_file)
+    # upload(s3, folder, output_vector)
 
+    name, ext = os.path.splitext(file)
+    output_file = f"{name}_processed{ext}"
+
+    output_cog_file = f"{name}_processed_COG{ext}"
+
+    output_vector_file = f"{name}_processed.parquet"
+    print(output_file)
+    print(output_cog_file)
+    print(output_vector_file)
+    # Reproject raster
+    if not exists_on_s3(s3, folder, output_file):
+        output_file = reproject_raster(file)
+        upload(s3, folder, output_file)
+    else:
+        print(f"{output_file} already exists on S3, skipping reprojection/upload.")
+
+    # Make COG
+    if not exists_on_s3(s3, folder, output_cog_file):
+        output_cog_file = make_cog(output_file)
+        upload(s3, folder, output_cog_file)
+    else:
+        print(f"{output_cog_file} already exists on S3, skipping COG conversion/upload.")
+
+    # Vectorize raster
+    if not exists_on_s3(s3, folder, output_vector_file):
+        output_vector_file, cols = make_vector(output_file)
+        upload(s3, folder, output_vector_file)
+    else:
+        print(f"{output_vector_file} already exists on S3, skipping vectorization/upload.")
+        # We still need column names
+        gdf = gpd.read_parquet(output_vector_file)
+        cols = gdf.drop('geom', axis=1).columns.to_list()
+    return cols
+    
 def reproject_raster(input_file, crs="EPSG:3310"):
     """
     Reproject rasters 
@@ -147,7 +184,7 @@ def make_vector(input_file, crs="EPSG:4326"):
 
     gdf.to_parquet(output_file)
     print(gdf)
-    return output_file
+    return output_file, gdf.drop('geom',axis = 1).columns.to_list()
 
 def filter_raster(s3, folder, file, percentile):
     """
@@ -168,31 +205,33 @@ def filter_raster(s3, folder, file, percentile):
     profile.update(dtype=rasterio.float64)
     with rasterio.open(new_file, "w", **profile) as dst:
         dst.write(filtered, 1)
-    process_raster(s3, folder, file)
-    return
+    cols = process_raster(s3, folder, file)
+    return cols
 
-def process_raster(s3, folder, file, file_name = None):
-    """
-    Driver function to process rasters 
-    """
-    if file_name:
-        file = file_name
-    output_file = reproject_raster(file)
-    upload(s3, folder, output_file)
-    output_cog_file = make_cog(output_file)
-    upload(s3, folder, output_cog_file)
-    output_vector = make_vector(output_file)
-    upload(s3, folder, output_vector)
-    return
-    
-def convert_pmtiles(folder, file):
+
+def convert_pmtiles(con, s3, folder, file):
     """
     Convert to PMTiles with tippecanoe 
     """
     name, ext = os.path.splitext(file)
     if ext != '.geojson':
-        con.read_parquet(file).execute().set_crs('epsg:3310').to_crs('epsg:4326').to_file(name+'.geojson')
+            (con.read_parquet(file).execute().set_crs('epsg:3310')
+             .to_crs('epsg:4326').to_file(name+'.geojson'))
     to_pmtiles(name+'.geojson', name+'.pmtiles', options = ['--extend-zooms-if-still-dropping'])
     upload(s3, folder, name+'.pmtiles')
     return
 
+def exists_on_s3(s3, folder, file):
+    """
+    Check if a file exists on S3
+    """
+    bucket, path = info(folder, file)
+    print(bucket)
+    print(path)
+    
+    try:
+        s3.stat_object(bucket, path)
+        return True
+    except S3Error:
+        return False
+