removing h3 code, will revisit later

later: fix h3 to raster, process h3 by looping.

preprocess/CBN-data.ipynb

@@ -18,11 +18,11 @@
     "from cng.utils import *\n",
     "from utils import *\n",
     "from h3_utils import *\n",
+    "s3 = s3_client()\n",
     "\n",
     "import os\n",
     "os.chdir('../data/')\n",
     "\n",
-    "s3 = s3_client()\n",
     "duckdb_install_h3()"
    ]
   },
@@ -48,11 +48,11 @@
     "folder = 'Counties'\n",
     "name = 'CA_counties'\n",
     "\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",
+    "unzip(s3, folder = folder, file = '30x30_Counties.zip')\n",
+    "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 = cols, zoom = 8)"
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -76,10 +76,11 @@
     "\n",
     "folder = 'Climate_zones'\n",
     "name = 'climate_zones_10'\n",
-    "# download(s3, 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,\n",
-    "          zoom = 8)\n"
+    "download(s3, folder = folder, file = 'clusters_10.tif')\n",
+    "process_raster(s3, folder = folder, file = 'clusters_10.tif', file_name = f\"{name}.tif\")\n",
+    "\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols,\n",
+    "          # zoom = 8)"
    ]
   },
   {
@@ -104,10 +105,10 @@
     "folder = 'Ecoregion'\n",
     "name = 'ACE_ecoregions'\n",
     "\n",
-    "# unzip(s3, folder = folder, file = '30x30_Ecoregions.zip')\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\")\n",
+    "unzip(s3, folder = folder, file = '30x30_Ecoregions.zip')\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.shp\")\n",
     "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -139,21 +140,21 @@
     "\n",
     "folder = 'Habitat'\n",
     "name = 'fveg22_1'\n",
-    "# unzip(s3, 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",
-    "cols = process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
-    "# upload(folder = folder, file = f'{name}_processed.tif.aux.xml')\n",
+    "subprocess.run(command, check=True)\n",
+    "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 = cols,\n",
-    "          zoom = 8)"
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols,\n",
+    "#           zoom = 8)"
    ]
   },
   {
@@ -186,14 +187,15 @@
     "folder = 'ACE_biodiversity'\n",
     "name = 'ACE_terrestrial_biodiversity_summary_ds2739'\n",
     "\n",
-    "# download(s3, folder = folder, file = 'Terrestrial_Biodiversity_Summary_-_ACE_[ds2739].geojson',\n",
-    "         # file_name = f\"{name}.geojson\")\n",
+    "download(s3, folder = folder, file = 'Terrestrial_Biodiversity_Summary_-_ACE_[ds2739].geojson',\n",
+    "         file_name = f\"{name}.geojson\")\n",
     "\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.geojson\")\n",
-    "cols = [item for item in cols if item not in [\"Hex_ID\",\"Shape__Area\",\"Shape__Length\"]]\n",
-    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.geojson\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)\n",
-    "# gdf = gpd.read_parquet(f\"{name}.parquet\")\n"
+    "process_vector(s3, folder = folder, file = f\"{name}.geojson\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.geojson\")\n",
+    "gdf = gpd.read_parquet(f\"{name}.parquet\")\n",
+    "\n",
+    "# cols = [item for item in cols if item not in [\"Hex_ID\",\"Shape__Area\",\"Shape__Length\"]]\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -222,7 +224,7 @@
     "            '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",
-    "    cols = process_vector(s3, folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = rank_df)\n",
+    "    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"
    ]
   },
@@ -261,7 +263,7 @@
     "        percentile = 0.95\n",
     "    threshold = gdf[col].quantile(percentile)\n",
     "    ace = gdf[gdf[col]>=threshold][cols]\n",
-    "    cols = process_vector(s3, folder = 'ACE_biodiversity/'+name, file = name+'.parquet',gdf = ace)\n",
+    "    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",
@@ -299,9 +301,10 @@
     "folder = 'Biodiversity_unique/Plant_richness'\n",
     "name = 'species_D'\n",
     "\n",
-    "# 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, zoom = 8)"
+    "download(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "filter_raster(s3, folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
+    "\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -327,8 +330,9 @@
     "name = 'endemicspecies_E'\n",
     "\n",
     "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, zoom = 8)"
+    "filter_raster(s3, folder = folder, file = f\"{name}.tif\", percentile = 80)\n",
+    "\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -370,9 +374,10 @@
     "folder = 'Connectivity_resilience/Resilient_connected_network_allcategories'\n",
     "name = 'rcn_wIntactBioCat_caOnly_2020-10-27'\n",
     "\n",
-    "cols = process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols, \n",
-    "           zoom = 8)"
+    "process_raster(s3, folder = folder, file = f\"{name}.tif\")\n",
+    "\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}_processed.parquet\", cols = cols, \n",
+    "           # zoom = 8)"
    ]
   },
   {
@@ -459,10 +464,11 @@
     "wetlands = ['Freshwater Emergent Wetland', 'Freshwater Forested/Shrub Wetland', 'Estuarine and Marine Wetland']\n",
     "gdf = gdf[gdf['WETLAND_TYPE'].isin(wetlands)]\n",
     "\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
-    "cols = [item for item in cols if item not in ['ACRES','Shape_Length','Shape_Area','__index_level_0__']]\n",
+    "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, zoom = 8)\n"
+    "\n",
+    "# cols = [item for item in cols if item not in ['ACRES','Shape_Length','Shape_Area','__index_level_0__']]\n",
+    "# geom_to_h3(con, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -566,28 +572,29 @@
     "\n",
     "folder = 'NBS_agriculture/Farmland_all'\n",
     "name = 'Important_Farmland_2018'\n",
-    "# unzip(s3, folder = folder, file = f\"{name}.zip\")\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.gdb\",crs = \"epsg:4326\")\n",
-    "cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "unzip(s3, folder = folder, file = f\"{name}.zip\")\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.gdb\",crs = \"epsg:4326\")\n",
+    "\n",
+    "convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n",
     "\n",
-    "# convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)\n",
+    "# cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)\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",
-    "# 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",
+    "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(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",
-    "# 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"
+    "gdf_grazing = gdf[gdf['polygon_ty'] == 'G']\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf_grazing)\n",
+    "convert_pmtiles(con, s3, folder = folder, file =f\"{name}.parquet\")\n"
    ]
   },
   {
@@ -640,15 +647,14 @@
     "\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_']>=2014]\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
-    "cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "gdf = gdf[~gdf['YEAR_'].isna()]\n",
+    "gdf['YEAR_'] = gdf['YEAR_'].astype('int64')\n",
+    "gdf = gdf[gdf['YEAR_']>=2014]\n",
+    "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_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
-    "\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
+    "# cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -732,13 +738,12 @@
     "folder = 'Progress_data_new_protection/Newly_counted_lands'\n",
     "name = 'newly_counted_lands_2024'\n",
     "\n",
+    "unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.shp\",crs = \"epsg:4326\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
-    "# unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.shp\",crs = \"epsg:4326\")\n",
-    "cols = [item for item in cols if item not in ['Shape_Leng', 'Shape_Area']]\n",
-    "\n",
-    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)\n"
+    "# cols = [item for item in cols if item not in ['Shape_Leng', 'Shape_Area']]\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
   {
@@ -764,8 +769,8 @@
     "name = 'DAC_2022'\n",
     "\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"
+    "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\")"
    ]
   },
   {
@@ -795,14 +800,22 @@
     "       .mutate(id=ibis.row_number().over()) #making a unique id \n",
     "      ).execute().set_crs('EPSG:3857')\n",
     "\n",
-    "cols = process_vector(s3, folder = folder, file = 'Priority Populations 4.0 Combined Layer.gdb',\n",
+    "process_vector(s3, folder = folder, file = 'Priority Populations 4.0 Combined Layer.gdb',\n",
     "               file_name = f\"{name}.parquet\", gdf = gdf)\n",
-    "cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
-    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)\n"
+    "# cols = [item for item in cols if item not in ['Shape_Length','Shape_Area']]\n",
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 8)"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "df1a939c-cb89-4a2f-8309-2819fe52ac45",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "markdown",
    "id": "a919ff5f-dff3-4db7-81c2-694f07f37d1d",
@@ -828,7 +841,7 @@
     "\n",
     "gdf = gpd.read_file('Priority Populations 4.0 Combined Layer.gdb')\n",
     "gdf = gdf[gdf['Designatio'] =='Low-income community']\n",
-    "cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
     "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")"
    ]
   },
@@ -853,42 +866,12 @@
     "folder = 'Progress_data_new_protection/Land_Status_Zone_Ecoregion_Counties'\n",
     "name = 'all_regions_reGAP_county_eco'\n",
     "\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",
+    "unzip(s3, folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp.zip')\n",
+    "process_vector(s3, folder = folder, file = 'Land_Status_Zone_Ecoregion_Counties.shp',\n",
     "               file_name = f\"{name}.parquet\")\n",
-    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 5)\n",
-    "\n",
-    "# convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "df6e2e1e-b74f-4b14-8140-7e425a3dec20",
-   "metadata": {},
-   "source": [
-    "# CA Nature data"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "16f9f330-c10c-4cec-9eba-0878aab9a5f7",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "%%time \n",
-    "con = ibis.duckdb.connect('ca_30x30_base',extensions = [\"spatial\", \"h3\"])\n",
-    "set_secrets(con)\n",
-    "\n",
-    "# file = 'ca-30x30-base.parquet'\n",
-    "folder = \"CA_Nature/2024/Preprocessing\"\n",
-    "name = 'ca-30x30-base'\n",
-    "# download(s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
     "\n",
-    "# cols = process_vector(s3, folder = folder, file = f\"{name}.parquet\", crs=\"EPSG:4326\")\n",
-    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, \n",
-    "           # zoom = 8)\n",
-    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n"
+    "# convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols = cols, zoom = 5)"
    ]
   },
   {
@@ -912,18 +895,18 @@
     "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, zoom = 8)\n",
+    "unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:3310\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "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, zoom = 8)\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, zoom = 8)"
+    "unzip(s3, folder = folder, file = f\"{name}.shp.zip\")\n",
+    "process_vector(s3, folder = folder, file = f\"{name}.shp\", crs=\"EPSG:3310\")\n",
+    "convert_pmtiles(con, s3, folder = folder, file = f\"{name}.parquet\")\n",
+    "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, zoom = 8)"
    ]
   }
  ],

preprocess/h3_utils.py

@@ -1,14 +1,13 @@
 from utils import *
 import re
 
-def convert_h3(con, s3, folder, file, cols, zoom):
+def convert_h3(con, s3, folder, file, cols, zoom, base_folder = "CBN/"):
     """
     Driver function to convert geometries to h3.
     If no zoom levels exist -> compute from geometry at target zoom.
-    If lower zoom exists -> compute children from max available until target zoom.
     """
     cols = ", ".join(cols) if isinstance(cols, list) else cols
-    bucket, path = info(folder, file)
+    bucket, path = info(folder, file, base_folder)
     path, file = os.path.split(path)
     name, ext = os.path.splitext(file)
     name = name.replace('-', '')
@@ -23,7 +22,7 @@ def convert_h3(con, s3, folder, file, cols, zoom):
             zooms.append(int(match.group(1)))
             
     if not zooms: # if no h3 files exist
-        print(f'No h3 files exists, computing {zoom} from geometry.')
+        print(f'No h3 files exists, computing zoom level {zoom} from geometry.')
         con.read_parquet(f"s3://{bucket}/{path}/{file}", table_name=name)
         h3_from_geom(con, name, cols, zoom)
         con.sql(f'''
@@ -38,21 +37,21 @@ def convert_h3(con, s3, folder, file, cols, zoom):
             print(f'Zoom {zoom} already exists!')
             return 
             
-        elif current_zoom < zoom: #compute child of most refined zoom level
-            print(f'Reading zoom {current_zoom}')
-            con.read_parquet(
-                f"s3://{bucket}/{path}/hex/zoom{current_zoom}/{name}.parquet",
-                table_name=f"h3_h{current_zoom}"
-            )
-            print(f'Computing {zoom} from {current_zoom}')
+        # elif current_zoom < zoom: #compute child of most refined zoom level
+        #     print(f'Reading zoom {current_zoom}')
+        #     con.read_parquet(
+        #         f"s3://{bucket}/{path}/hex/zoom{current_zoom}/{name}.parquet",
+        #         table_name=f"h3_h{current_zoom}"
+        #     )
+        #     print(f'Computing {zoom} from {current_zoom}')
             
-            for z in range(current_zoom + 1, zoom + 1):
-                print(f'Current zoom {z}')
-                h3_from_parent(con, z)
-                con.sql(f'''
-                    SELECT *, UNNEST(h3_cell_to_children(h{z-1}, {z})) AS h{z}
-                    FROM h3_h{z-1}
-                ''').to_parquet(f"s3://{bucket}/{path}/hex/zoom{z}/{name}.parquet")
+        #     for z in range(current_zoom + 1, zoom + 1):
+        #         print(f'Current zoom {z}')
+        #         h3_from_parent(con, z)
+        #         con.sql(f'''
+        #             SELECT *, UNNEST(h3_cell_to_children(h{z-1}, {z})) AS h{z}
+        #             FROM h3_h{z-1}
+        #         ''').to_parquet(f"s3://{bucket}/{path}/hex/zoom{z}/{name}.parquet")
   
     
 def h3_from_geom(con, name, cols, zoom):
@@ -60,20 +59,19 @@ def h3_from_geom(con, name, cols, zoom):
     Computes hexes directly from geometry.
     """
     con.raw_sql(f'''
-        CREATE OR REPLACE TEMP TABLE t2 AS
-        WITH t1 AS (
-            SELECT {cols}, ST_Dump(geom) AS geom 
-            FROM {name}
-        )
-        SELECT {cols},
-               h3_polygon_wkt_to_cells_string(ST_Force2D(UNNEST(geom).geom), {zoom}) AS h{zoom}
-        FROM t1
+    CREATE OR REPLACE TEMP TABLE t2 AS
+    SELECT {cols},
+           h3_polygon_wkt_to_cells_string(ST_Force2D(dump.geom), {zoom}) AS h{zoom}
+    FROM (
+        SELECT {cols}, UNNEST(ST_Dump(geom)) AS dump
+        FROM {name}
+    )
     ''')
 
 
-def h3_from_parent(con, zoom):
-    con.raw_sql(f'''
-        CREATE OR REPLACE TEMP TABLE h3_h{zoom} AS
-        SELECT *, UNNEST(h3_cell_to_children(h{zoom-1}, {zoom})) AS h{zoom}
-        FROM h3_h{zoom-1}
-    ''')
+# def h3_from_parent(con, zoom):
+#     con.raw_sql(f'''
+#         CREATE OR REPLACE TEMP TABLE h3_h{zoom} AS
+#         SELECT *, UNNEST(h3_cell_to_children(h{zoom-1}, {zoom})) AS h{zoom}
+#         FROM h3_h{zoom-1}
+#     ''')

preprocess/utils.py

@@ -10,48 +10,50 @@ import ibis
 from ibis import _
 
 import rasterio
-from rasterio.features import shapes
-from shapely.geometry import shape
+from rasterio.transform import xy
+from shapely.geometry import Point
 import numpy as np
+from pyproj import Transformer
 
-
-def info(folder, file, bucket = "public-ca30x30", base_folder = 'CBN/'):
+def info(folder, file, base_folder, bucket = "public-ca30x30"):
     """
     Extract minio path to upload/download data 
     """
-    path = os.path.join(base_folder, folder, file)
+    if (folder is None) & (base_folder is None):
+        path = file
+    else:
+        path = os.path.join(base_folder, folder, file)
     # path = os.path.join(folder, file)
     return bucket, path 
     
-def download(s3, folder, file, file_name = None):
+def download(s3, folder, file, file_name = None, base_folder = "CBN/"):
     """
     Downloading file from minio
     """
     if not file_name: 
         file_name = file
-    bucket, path = info(folder, file)
-    s3.fget_object(bucket, path ,file_name) 
+    bucket, path = info(folder, file, base_folder)
+    s3.fget_object(bucket, path , file_name) 
     return
 
-def upload(s3, folder, file):
+def upload(s3, folder, file, base_folder = "CBN/"):
     """
     Uploading file from minio
     """
-    bucket, path = info(folder, file)
+    bucket, path = info(folder, file, base_folder)
     s3.fput_object(bucket, path ,file) 
     return
 
-def unzip(s3, folder, file):
+def unzip(s3, folder, file, base_folder = "CBN/"):
     """
     Unzipping zip files 
     """
-    download(s3, folder, file)
+    download(s3, folder, file, base_folder)
     with zipfile.ZipFile(file, 'r') as zip_ref:
         zip_ref.extractall()
     return 
 
-# 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"):
+def process_vector(s3, folder, file, file_name = None, gdf = None, crs="EPSG:4326", base_folder = "CBN/"):
     """
     Driver function to process vectors 
     """
@@ -67,11 +69,11 @@ def process_vector(s3, folder, file, file_name = None, gdf = None, crs="EPSG:432
     name, ext = os.path.splitext(file)
     parquet_file = f"{name}{'.parquet'}"
     gdf.to_parquet(parquet_file)
-    upload(s3, folder, parquet_file)
-
-    return gdf.drop('geom',axis = 1).columns.to_list()
+    upload(s3, folder, parquet_file, base_folder)
+    # return gdf.drop('geom',axis = 1).columns.to_list()
+    return
 
-def process_raster(s3, folder, file, file_name = None):
+def process_raster(s3, folder, file, file_name = None, base_folder = "CBN/"):
     """
     Driver function to process rasters 
     """
@@ -82,29 +84,30 @@ def process_raster(s3, folder, file, file_name = None):
     output_cog_file = f"{name}_processed_COG{ext}"
     output_vector_file = f"{name}_processed.parquet"
     # Reproject raster
-    if not exists_on_s3(s3, folder, output_file):
+    if not exists_on_s3(s3, folder, output_file, base_folder):
         output_file = reproject_raster(file)
-        upload(s3, folder, output_file)
+        upload(s3, folder, output_file, base_folder)
     else:
         print(f"{output_file} already exists on S3, skipping reprojection/upload.")
 
     # Make COG
-    if not exists_on_s3(s3, folder, output_cog_file):
+    if not exists_on_s3(s3, folder, output_cog_file, base_folder):
         output_cog_file = make_cog(output_file)
-        upload(s3, folder, output_cog_file)
+        upload(s3, folder, output_cog_file, base_folder)
     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
+    # # Vectorize raster
+    # if not exists_on_s3(s3, folder, output_vector_file, base_folder):
+    #     output_vector_file, cols = make_vector(output_file)
+    #     upload(s3, folder, output_vector_file, base_folder)
+    # 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
+    return 
     
 def reproject_raster(input_file, crs="EPSG:3310"):
     """
@@ -147,31 +150,37 @@ def make_cog(input_file, crs="EPSG:4326"):
         print(f"Error occurred during processing: {e}")
     return output_file 
 
-def make_vector(input_file, crs="EPSG:4326"):
-    """
-    Converting rasters to vector formats in order to convert to h3
-    """
-    name, ext = os.path.splitext(input_file)
-    output_file = f"{name}.parquet"
-    # Open raster
-    with rasterio.open(input_file) as src:
-        image = src.read(1)  # read first band
-        mask = image != src.nodata  # mask out nodata
+# def make_vector(input_file, crs="EPSG:4326"):
+#     """
+#     Converting rasters to vector formats in order to convert to h3
+#     """
+#     name, ext = os.path.splitext(input_file)
+#     output_file = f"{name}.parquet"
     
-        results = (
-            {"geom": shape(geom), "value": value}
-            for geom, value in shapes(image, mask=mask, transform=src.transform)
-        )
-    
-    gdf = gpd.GeoDataFrame.from_records(results)
-    gdf.set_geometry('geom', inplace=True)
-    gdf['id'] = np.arange(len(gdf))
-    gdf.set_crs(src.crs, inplace=True)
-    if gdf.crs != crs:
-        gdf.to_crs(crs, inplace=True)
-
-    gdf.to_parquet(output_file)
-    return output_file, gdf.drop('geom',axis = 1).columns.to_list()
+#     with rasterio.open(input_file) as src:
+#         band = src.read(1)  # read first band
+#         mask = band != src.nodata  # mask out nodata
+#         rows, cols = np.where(mask)
+#         x, y = rasterio.transform.xy(src.transform, rows, cols, offset = "center")
+
+#         # reproject 
+#         if src.crs and src.crs.to_string() != crs:
+#             transformer = Transformer.from_crs(src.crs, crs, always_xy=True)
+#             x, y = transformer.transform(x, y)
+#             crs_out = crs
+#         else:
+#             crs_out = src.crs
+
+#         gdf = gpd.GeoDataFrame(
+#             {"value": band[rows, cols]},
+#             geometry=[Point(xy) for xy in zip(x, y)],
+#             crs=crs_out
+#         )
+
+#     gdf.rename_geometry('geom', inplace=True)
+#     gdf['id'] = np.arange(len(gdf))
+#     gdf.to_parquet(output_file)
+#     return output_file, gdf.drop('geom',axis = 1).columns.to_list()
 
 def filter_raster(s3, folder, file, percentile):
     """
@@ -192,11 +201,11 @@ 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)
-    cols = process_raster(s3, folder, file)
-    return cols
-
+    process_raster(s3, folder, file)
+    # return cols
+    return
 
-def convert_pmtiles(con, s3, folder, file):
+def convert_pmtiles(con, s3, folder, file, base_folder = "CBN/"):
     """
     Convert to PMTiles with tippecanoe 
     """
@@ -205,14 +214,14 @@ def convert_pmtiles(con, s3, folder, file):
             (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')
+    upload(s3, folder, name+'.pmtiles', base_folder)
     return
 
-def exists_on_s3(s3, folder, file):
+def exists_on_s3(s3, folder, file, base_folder = "CBN/"):
     """
     Check if a file exists on S3
     """
-    bucket, path = info(folder, file)
+    bucket, path = info(folder, file, base_folder)
     try:
         s3.stat_object(bucket, path)
         return True