wip; h3!!!

preprocess/CBN-data.ipynb

@@ -15,9 +15,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from cng.utils import set_secrets, s3_client, to_pmtiles\n",
+    "from cng.utils import *\n",
+    "from utils import *\n",
+    "from h3_utils import *\n",
     "s3 = s3_client()\n",
     "\n",
+    "duckdb_install_h3()\n",
+    "from minio.error import S3Error\n",
+    "\n",
+    "\n",
     "import zipfile\n",
     "import os\n",
     "import subprocess\n",
@@ -26,140 +32,15 @@
     "import geopandas as gpd\n",
     "import ibis\n",
     "from ibis import _\n",
-    "con = ibis.duckdb.connect(extensions=[\"spatial\"])\n",
+    "# con = ibis.duckdb.connect('temp',extensions = [\"spatial\", \"h3\"])\n",
+    "# set_secrets(con)\n",
     "\n",
     "import rasterio\n",
+    "from rasterio.features import shapes\n",
+    "from shapely.geometry import shape\n",
     "import numpy as np"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "a6e4eac1-9bde-4d83-a386-60a055dfe402",
-   "metadata": {},
-   "source": [
-    "#### Helper functions"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "id": "0fc80b2e-375a-4295-95cf-5bfc313b78fd",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def info(folder, file, bucket = \"public-ca30x30\", base_folder = 'CBN-data/'):\n",
-    "    path = os.path.join(base_folder, folder, file)\n",
-    "    return bucket, path \n",
-    "    \n",
-    "def download(folder, file, file_name = None):\n",
-    "    if not file_name: \n",
-    "        file_name = file\n",
-    "    bucket, path = info(folder, file)\n",
-    "    s3.fget_object(bucket, path ,file_name) \n",
-    "    return\n",
-    "\n",
-    "def upload(folder, file):\n",
-    "    bucket, path = info(folder, file)\n",
-    "    s3.fput_object(bucket, path ,file) \n",
-    "    return\n",
-    "\n",
-    "def unzip(folder, file):\n",
-    "    download(folder, file)\n",
-    "    with zipfile.ZipFile(file, 'r') as zip_ref:\n",
-    "        zip_ref.extractall()\n",
-    "    return \n",
-    "\n",
-    "def upload_parquet(folder, file, gdf):\n",
-    "    name, ext = os.path.splitext(file)\n",
-    "    parquet_file = f\"{name}{'.parquet'}\"\n",
-    "    gdf.to_parquet(parquet_file)\n",
-    "    upload(folder, parquet_file)\n",
-    "    return  \n",
-    "\n",
-    "def process_vector(folder, file, file_name = None, gdf = None, crs=\"EPSG:3310\"):\n",
-    "    if gdf is None:\n",
-    "        gdf = gpd.read_file(file)\n",
-    "    gdf = gdf.to_crs(crs)\n",
-    "    gdf = gdf.rename_geometry('geom')\n",
-    "    if file_name:\n",
-    "        file = file_name\n",
-    "    upload_parquet(folder, file, gdf)\n",
-    "    return \n",
-    "\n",
-    "def reproject_raster(input_file, crs=\"EPSG:3310\"):\n",
-    "    suffix = '_processed'\n",
-    "    name, ext = os.path.splitext(input_file)\n",
-    "    output_file = f\"{name}{suffix}{ext}\"\n",
-    "    command = [\n",
-    "        \"gdalwarp\",\n",
-    "        \"-t_srs\", crs,\n",
-    "        input_file,\n",
-    "        output_file \n",
-    "        ]\n",
-    "    try:\n",
-    "        subprocess.run(command, check=True)\n",
-    "        print(f\"Reprojection successful!\")\n",
-    "    except subprocess.CalledProcessError as e:\n",
-    "        print(f\"Error occurred during reprojection: {e}\")\n",
-    "    return output_file \n",
-    "\n",
-    "def make_cog(input_file, crs=\"EPSG:4326\"):\n",
-    "    suffix = '_COG'\n",
-    "    name, ext = os.path.splitext(input_file)\n",
-    "    output_file = f\"{name}{suffix}{ext}\"\n",
-    "    command = [\n",
-    "        \"gdalwarp\",\n",
-    "        \"-t_srs\", crs,\n",
-    "        \"-of\", \"COG\",\n",
-    "        input_file,\n",
-    "        output_file \n",
-    "        ]\n",
-    "    try:\n",
-    "        subprocess.run(command, check=True)\n",
-    "        print(f\"Successful!\")\n",
-    "    except subprocess.CalledProcessError as e:\n",
-    "        print(f\"Error occurred during processing: {e}\")\n",
-    "    return output_file \n",
-    "\n",
-    "def process_raster(folder, file, file_name = None):\n",
-    "    if file_name:\n",
-    "        file = file_name\n",
-    "    output_file = reproject_raster(file)\n",
-    "    upload(folder, output_file)\n",
-    "    output_cog_file = make_cog(output_file)\n",
-    "    upload(folder, output_cog_file)\n",
-    "    return\n",
-    "\n",
-    "def filter_raster(folder, file, percentile):\n",
-    "    with rasterio.open(file) as src:\n",
-    "        data = src.read(1)  # Read the first band\n",
-    "        profile = src.profile\n",
-    "    \n",
-    "    # mask no data values\n",
-    "    masked_data = np.ma.masked_equal(data, src.nodata)\n",
-    "\n",
-    "    # compute percentile/threshold \n",
-    "    p = np.percentile(masked_data.compressed(),percentile)\n",
-    "    filtered = np.where(data >= p, data, src.nodata)\n",
-    "    \n",
-    "    name, ext = os.path.splitext(file)\n",
-    "    new_file = f\"{name}{'_'}{percentile}{'percentile'}{ext}\"\n",
-    "\n",
-    "    profile.update(dtype=rasterio.float64)\n",
-    "    with rasterio.open(new_file, \"w\", **profile) as dst:\n",
-    "        dst.write(filtered, 1)\n",
-    "        \n",
-    "    process_raster(folder, file)\n",
-    "    return\n",
-    "\n",
-    "def convert_pmtiles(folder, parquet_file):\n",
-    "    name, ext = os.path.splitext(parquet_file)\n",
-    "    con.read_parquet(parquet_file).execute().set_crs('epsg:3310').to_crs('epsg:4326').to_file(name+'.geojson')\n",
-    "    to_pmtiles(name+'.geojson', name+'.pmtiles', options = ['--extend-zooms-if-still-dropping'])\n",
-    "    upload(folder, name+'.pmtiles')\n",
-    "    return"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "2219c881-1017-4def-9d0c-c83b5541b5d2",
@@ -175,8 +56,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Counties', file = '30x30_Counties.zip')\n",
-    "process_vector(folder = 'Counties', file = 'CA_counties.shp')"
+    "%%time \n",
+    "con = ibis.duckdb.connect('counties',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "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",
+    "\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= \"COUNTY_NAM\")"
    ]
   },
   {
@@ -194,8 +85,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "download(folder = 'Climate_zones', file = 'clusters_10.tif')\n",
-    "process_raster(folder = 'Climate_zones', file = 'clusters_10.tif', file_name =  'climate_zones_10.tif')"
+    "%%time \n",
+    "con = ibis.duckdb.connect('climate_zones',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\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\")"
    ]
   },
   {
@@ -213,8 +111,17 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Ecoregion', file = '30x30_Ecoregions.zip')\n",
-    "process_vector(folder = 'Ecoregion', file = 'ACE_ecoregions.shp')"
+    "%%time \n",
+    "con = ibis.duckdb.connect('ecoregion',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "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",
+    "\n",
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= 'CA_Ecoregi')"
    ]
   },
   {
@@ -241,7 +148,7 @@
    "outputs": [],
    "source": [
     "# download(folder = 'Habitat', file = 'CWHR13_2022.tif')\n",
-    "# process_raster(folder = 'Habitat', file = 'CWHR13_2022.tif')"
+    "# process_raster(s3, folder = 'Habitat', file = 'CWHR13_2022.tif')"
    ]
   },
   {
@@ -251,18 +158,27 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Habitat', file = 'fveg221gdb.zip')\n",
+    "%%time\n",
+    "con = ibis.duckdb.connect('habitat',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "folder = 'Habitat'\n",
+    "name = 'fveg22_1'\n",
     "\n",
-    "command = [\n",
-    "    \"gdalwarp\",\n",
-    "    \"-of\", \"GTiff\",\n",
-    "    'fveg22_1.gdb',\n",
-    "    'fveg22_1.tif' \n",
-    "    ]\n",
+    "# unzip(folder = folder, file = 'fveg221gdb.zip')\n",
     "\n",
-    "subprocess.run(command, check=True)\n",
-    "process_raster(folder = 'Habitat', file = 'fveg22_1.tif')\n",
-    "upload(folder = 'Habitat', file = 'fveg22_1_processed.tif.aux.xml')\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",
+    "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\")"
    ]
   },
   {
@@ -288,12 +204,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "download(folder = 'ACE_biodiversity', file = 'Terrestrial_Biodiversity_Summary_-_ACE_[ds2739].geojson',\n",
-    "         file_name = 'ACE_biodiversity_all_ds2739.geojson')\n",
-    "gdf = gpd.read_file('ACE_biodiversity_all_ds2739.geojson')\n",
+    "%%time \n",
+    "folder = 'ACE_biodiversity'\n",
+    "name = 'ACE_terrestrial_biodiversity_summary_ds2739'\n",
+    "\n",
+    "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",
     "\n",
-    "to_pmtiles('ACE_biodiversity_all_ds2739.geojson', 'ACE_biodiversity_ds2739_all.pmtiles', options = ['--extend-zooms-if-still-dropping'])\n",
-    "upload(folder = 'ACE_biodiversity',file = 'ACE_biodiversity_ds2739_all.pmtiles')"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= \"OBJECTID\")\n",
+    "# gdf = gpd.read_parquet(f\"{name}.parquet\")\n"
    ]
   },
   {
@@ -323,7 +245,7 @@
     "    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, parquet_file = name+'.parquet')\n"
+    "    convert_pmtiles(folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n"
    ]
   },
   {
@@ -362,7 +284,7 @@
     "    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, parquet_file = name+'.parquet')\n",
+    "    convert_pmtiles(folder ='ACE_biodiversity/'+name, file = name+'.parquet')\n",
     "\n",
     "\n",
     "# calculate 80% percentile, filter to those >= threshold. \n",
@@ -392,8 +314,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "download(folder = 'Biodiversity_unique/Plant_richness', file = 'species_D.tif')\n",
-    "filter_raster(folder = 'Biodiversity_unique/Plant_richness', file = 'species_D.tif', percentile = 80)\n"
+    "%%time \n",
+    "con = ibis.duckdb.connect('plant',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "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\")"
    ]
   },
   {
@@ -411,10 +341,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "download(folder = 'Biodiversity_unique/Rarityweighted_endemic_plant_richness', file = 'endemicspecies_E.tif')\n",
+    "%%time \n",
+    "con = ibis.duckdb.connect('end_plant',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "folder = 'Biodiversity_unique/Rarityweighted_endemic_plant_richness'\n",
+    "name = 'endemicspecies_E'\n",
     "\n",
-    "filter_raster(folder = 'Biodiversity_unique/Rarityweighted_endemic_plant_richness',\n",
-    "               file = 'endemicspecies_E.tif', percentile = 80)"
+    "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\")"
    ]
   },
   {
@@ -449,8 +385,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "process_raster(folder = 'Connectivity_resilience/Resilient_connected_network_allcategories',\n",
-    "               file = 'rcn_wIntactBioCat_caOnly_2020-10-27.tif')"
+    "%%time \n",
+    "con = ibis.duckdb.connect('CRN',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "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\")"
    ]
   },
   {
@@ -524,14 +467,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Freshwater_resources/Wetlands', file = 'CA_geodatabase_wetlands.zip')\n",
+    "%%time \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",
     "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",
-    "process_vector(folder = 'Freshwater_resources/Wetlands', file = 'CA_wetlands.parquet', gdf = gdf)\n",
-    "convert_pmtiles(folder ='Freshwater_resources/Wetlands', parquet_file ='CA_wetlands.parquet')\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"
    ]
   },
   {
@@ -625,23 +575,39 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "7e39ee64-3123-43f0-9552-319010708b34",
+   "id": "fdf7061f-7598-4303-bb77-38f836feac8a",
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'NBS_agriculture/Farmland', file = 'Important_Farmland_2018.zip')\n",
+    "%%time \n",
+    "\n",
+    "folder = 'NBS_agriculture/Farmland'\n",
+    "unzip(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",
     "\n",
     "# only pick a subset \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 = 'NBS_agriculture/Farmland', file = 'Farmland_2018.parquet', gdf = gdf_farmland)\n",
-    "convert_pmtiles(folder ='NBS_agriculture/Farmland', parquet_file ='Farmland_2018.parquet')\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_grazing = gdf[gdf['polygon_ty'] == 'G']\n",
-    "process_vector(folder = 'NBS_agriculture/Lands_suitable_grazing', \n",
-    "               file = 'Grazing_land_2018.parquet', gdf = gdf_grazing)\n",
-    "convert_pmtiles(folder ='NBS_agriculture/Lands_suitable_grazing', parquet_file ='Grazing_land_2018.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"
    ]
   },
   {
@@ -678,6 +644,14 @@
     "Only YEAR >= 2014. "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "425f9149-d8ac-437a-9572-301bd1b1bec8",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -685,14 +659,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Climate_risks/Historical_fire_perimeters', file = 'fire23-1gdb.zip')\n",
+    "%%time \n",
+    "folder = 'Climate_risks/Historical_fire_perimeters'\n",
+    "name = 'calfire_2023'\n",
+    "\n",
+    "unzip(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",
+    "# 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",
     "\n",
-    "process_vector(folder = 'Climate_risks/Historical_fire_perimeters', file = 'calfire_2023.parquet', gdf = gdf)\n",
-    "convert_pmtiles(folder ='Climate_risks/Historical_fire_perimeters', parquet_file ='calfire_2023.parquet')\n"
+    "convert_h3(con, s3, folder = folder, file = f\"{name}.parquet\", cols= ['INC_NUM','FIRE_NAME','YEAR_'])"
    ]
   },
   {
@@ -769,9 +748,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Progress_data_new_protection/Newly_counted_lands', file = 'newly_counted_lands_2024.shp.zip')\n",
-    "process_vector(folder = 'Progress_data_new_protection/Newly_counted_lands', file = 'newly_counted_lands_2024.shp')\n",
-    "convert_pmtiles(folder ='Progress_data_new_protection/Newly_counted_lands', parquet_file ='newly_counted_lands_2024.parquet')\n"
+    "%%time \n",
+    "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"
    ]
   },
   {
@@ -789,10 +773,14 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Progress_data_new_protection/DAC', file = 'sb535dacgdbf2022gdb.zip')\n",
-    "process_vector(folder = 'Progress_data_new_protection/DAC', file = 'SB535DACgdb_F_2022.gdb',\n",
-    "               file_name = 'DAC_2022.parquet')\n",
-    "convert_pmtiles(folder ='Progress_data_new_protection/DAC', parquet_file ='DAC_2022.parquet')\n"
+    "%%time \n",
+    "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"
    ]
   },
   {
@@ -809,7 +797,42 @@
    "id": "fc236b85-9fc4-4589-8dd9-5efb7f2e9614",
    "metadata": {},
    "outputs": [],
-   "source": []
+   "source": [
+    "%%time \n",
+    "con = ibis.duckdb.connect('priority_pop',extensions = [\"spatial\", \"h3\"])\n",
+    "set_secrets(con)\n",
+    "\n",
+    "folder = 'Progress_data_new_protection/Priority_populations'\n",
+    "name = 'CalEnviroScreen4'\n",
+    "unzip(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",
+    "               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()"
+   ]
   },
   {
    "cell_type": "markdown",
@@ -826,16 +849,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Progress_data_new_protection/Low_income_communities',\n",
-    "      file = 'Priority Populations 4.0 Geodatabase.zip')\n",
+    "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",
     "\n",
     "gdf = gpd.read_file('Priority Populations 4.0 Combined Layer.gdb')\n",
     "gdf = gdf[gdf['Designatio'] =='Low-income community']\n",
-    "\n",
-    "\n",
-    "process_vector(folder = 'Progress_data_new_protection/Low_income_communities', \n",
-    "               file = 'low_income_CalEnviroScreen4.parquet',gdf = gdf)\n",
-    "convert_pmtiles(folder ='Progress_data_new_protection/Low_income_communities', parquet_file ='low_income_CalEnviroScreen4.parquet')\n"
+    "process_vector(folder = folder, file = f\"{name}.parquet\", gdf = gdf)\n",
+    "convert_pmtiles(folder = folder, file = f\"{name}.parquet\")"
    ]
   },
   {
@@ -853,12 +875,117 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "unzip(folder = 'Progress_data_new_protection/Land_Status_Zone_Ecoregion_Counties',\n",
-    "      file = 'Land_Status_Zone_Ecoregion_Counties.shp.zip')\n",
-    "process_vector(folder = 'Progress_data_new_protection/Land_Status_Zone_Ecoregion_Counties', file = 'Land_Status_Zone_Ecoregion_Counties.shp',\n",
-    "               file_name = 'all_regions_reGAP_county_eco.parquet')\n",
-    "convert_pmtiles(folder ='Progress_data_new_protection/Land_Status_Zone_Ecoregion_Counties', parquet_file ='all_regions_reGAP_county_eco.parquet')\n"
+    "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",
+    "               file_name = f\"{name}.parquet\")\n",
+    "convert_pmtiles(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": "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,
+   "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 = \"Preprocessing\"\n",
+    "name = 'ca-30x30-base'\n",
+    "# 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"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "908786ec-2a86-4fb0-a47a-9de364254806",
+   "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"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65d98aa3-041f-42d6-8448-6fa8a05f5850",
+   "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()"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa960a99-3c79-4e67-becf-a1cb397aa5fb",
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -877,7 +1004,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.8"
+   "version": "3.12.9"
   }
  },
  "nbformat": 4,

preprocess/h3_utils.py

@@ -0,0 +1,189 @@
+
+def geom_to_h3(con, name, cols, zoom):
+    """
+    Computes hexes    
+    """
+    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(UNNEST(geom).geom, {zoom}) AS h{zoom}
+        FROM t1
+    ''')
+    
+def check_size(con, name, zoom, sample_size=100):
+    """
+    Estimating size of geoms to decide if we need to process in chunks 
+    """
+    query = f"""
+        SELECT 
+            avg(len(h3_polygon_wkt_to_cells_string(ST_AsText(geom), {zoom}))::DOUBLE) AS avg_h3_len,
+            max(len(h3_polygon_wkt_to_cells_string(ST_AsText(geom), {zoom}))) AS max_h3_len,
+            count(*) AS total_rows
+        FROM {name}
+        USING SAMPLE {sample_size}
+    """
+    stats = con.sql(query).execute()
+    avg_len = stats.iloc[0]['avg_h3_len']
+    max_len = stats.iloc[0]['max_h3_len']
+    total_rows = con.table(name).count().execute()
+
+    est_total_h3 = avg_len * total_rows
+
+    print(f"Estimated total H3 cells: {est_total_h3:,.0f}")
+    print(f"Max H3 cells in one geometry: {max_len:,}")
+
+    return est_total_h3, max_len
+
+
+def write_large_geoms(con, s3, bucket, path, name, zoom="8", geom_len_threshold=10_000):
+    """
+    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}")
+
+        #  check if file already exists in minio
+        try:
+            s3.stat_object(bucket, large_key)
+            print(f"⏩ Skipping existing large batch: {large_key}")
+            offset += limit
+            i += 1
+            continue
+        except S3Error as err:
+            if err.code != "NoSuchKey":
+                raise  
+
+        print(f"📝 Writing large geometry batch {i} → {large_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}
+        ''')
+
+        q.to_parquet(f"s3://{bucket}/{large_key}")
+
+        if q.count().execute() == 0:
+            break
+
+        offset += limit
+        i += 1
+
+    return i
+
+    
+def join_large_geoms(con, s3, bucket, path, name):
+    """
+    If we had to process large geoms individually, join those datasets after conversion.
+    """
+    # 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
+    # join if it exists 
+    con.raw_sql(f'''
+        COPY (
+            SELECT * FROM read_parquet('s3://{bucket}/{path}/hex/{name}_large_*.parquet')
+        )
+        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):
+    """
+    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)
+            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'''
+            SELECT *, UNNEST(h{zoom}) AS h{zoom}
+            FROM t2
+            WHERE len(h{zoom}) <= {geom_len_threshold}
+            LIMIT {limit} OFFSET {offset}
+        ''')
+        q.to_parquet(f"s3://{bucket}/{chunk_path}")
+        if q.count().execute() == 0:
+            break
+        offset += limit
+        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)
+    join_large_geoms(con, s3, bucket, path, name)
+    return i
+
+
+
+def join_chunked(bucket, path, name):
+    """
+    If we had to chunk the data, join those datasets after conversion.
+    """
+    con.raw_sql(f'''
+        COPY (
+        SELECT * FROM read_parquet('s3://{bucket}/{path}/hex/{name}_chunk*.parquet')
+        )
+        TO 's3://{bucket}/{path}/hex/{name}.parquet'
+        (FORMAT PARQUET)
+        ''')
+
+# 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):
+    """
+    Driver function to convert geometries to h3
+    """
+    cols = ", ".join(cols) if isinstance(cols, list) else cols
+    bucket, path = info(folder, file)
+    path, file = os.path.split(path)
+    name, ext = os.path.splitext(file)
+    name = name.replace('-', '')
+
+    print(f"Processing: {name}")
+    con.read_parquet(f"s3://{bucket}/{path}/{file}", table_name=name)
+
+    # Decide to chunk or not
+    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:
+        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)
+        join_chunked(con, bucket, path, name)
+    else:
+        print("Writing single output")
+        geom_to_h3(con, name, cols, zoom)
+        con.sql(f'''
+            SELECT *, UNNEST(h{zoom}) AS h{zoom}
+            FROM t2
+        ''').to_parquet(f"s3://{bucket}/{path}/hex/{name}.parquet")

preprocess/preprocess.ipynb

@@ -10,7 +10,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "id": "f7e6298c-d886-432a-a1b7-c3fee914c24f",
    "metadata": {
     "editable": true,
@@ -25,12 +25,13 @@
     "data_path = '../data/CBN-layers/'\n",
     "os.chdir(data_path)\n",
     "\n",
-    "from cng.utils import ST_MakeValid, set_secrets, s3_client, s3_cp, to_pmtiles\n",
+    "from cng.utils import set_secrets, s3_client, s3_cp, to_pmtiles\n",
     "s3 = s3_client()\n",
     "\n",
+    "     \n",
     "import ibis\n",
-    "\n",
     "from ibis import _\n",
+    "import ibis.expr.datatypes as dt\n",
     "con = ibis.duckdb.connect(extensions=[\"spatial\"])\n",
     "\n",
     "import geopandas as gpd\n",
@@ -47,7 +48,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "id": "63dd33b8-6d3c-4852-9899-6ed5775d19c0",
    "metadata": {},
    "outputs": [],
@@ -60,7 +61,13 @@
     "    else:\n",
     "        path = os.path.join(bucket,base_folder,folder,file)\n",
     "    url = minio+path\n",
-    "    return url"
+    "    return url\n",
+    "\n",
+    "\n",
+    "# usage: t.mutate(geom_valid = ST_MakeValid(t.geom))\n",
+    "@ibis.udf.scalar.builtin\n",
+    "def ST_MakeValid(geom) -> dt.geometry:\n",
+    " ..."
    ]
   },
   {
@@ -73,14 +80,14 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 3,
    "id": "13214bbe-3a74-4247-981f-5a6eb6c486f5",
    "metadata": {},
    "outputs": [],
    "source": [
     "# CA Nature data \n",
-    "# ca_raw_parquet = \"https://data.source.coop/cboettig/ca30x30/ca_areas.parquet\"\n",
-    "ca_raw_parquet = 'ca_areas.parquet'\n",
+    "ca_raw_parquet = \"https://data.source.coop/cboettig/ca30x30/ca_areas.parquet\"\n",
+    "# ca_raw_parquet = 'ca_areas.parquet'\n",
     "\n",
     "# Boundary of CA, used to computed 'non-conserved' areas\n",
     "ca_boundary_parquet = get_url('Preprocessing','ca_boundary.parquet',base_folder = None)\n",
@@ -160,11 +167,45 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 7,
    "id": "0f9666d1-7c2b-45af-9399-e4189bba34f5",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "52ef18913f17417299860d91e36e9dbd",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "FloatProgress(value=0.0, layout=Layout(width='auto'), style=ProgressStyle(bar_color='black'))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "CPU times: user 4min 28s, sys: 6.1 s, total: 4min 34s\n",
+      "Wall time: 2min 18s\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "<minio.helpers.ObjectWriteResult at 0x7ff0943c7710>"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
+    "%%time \n",
     "# match CA Nature schema \n",
     "\n",
     "non_conserved = (con.read_parquet(ca_nonconserved_url)\n",
@@ -200,7 +241,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 4,
    "id": "a3d4f189-1563-4868-9f1f-64d67569df27",
    "metadata": {},
    "outputs": [],
@@ -257,7 +298,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 5,
    "id": "a59c976b-3c36-40f9-a15b-cefcd155c647",
    "metadata": {},
    "outputs": [],
@@ -303,11 +344,60 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 6,
    "id": "4d6177e2-8ece-4eb9-acc2-5fb5c5beb8bb",
    "metadata": {},
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "09f24f1359a84ae2a4b69360cc8e852b",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "FloatProgress(value=0.0, layout=Layout(width='auto'), style=ProgressStyle(bar_color='black'))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "data": {
+      "application/vnd.jupyter.widget-view+json": {
+       "model_id": "c10ce980d24e45b6bad9b8a70c176f2c",
+       "version_major": 2,
+       "version_minor": 0
+      },
+      "text/plain": [
+       "FloatProgress(value=0.0, layout=Layout(width='auto'), style=ProgressStyle(bar_color='black'))"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/opt/conda/lib/python3.12/site-packages/ibis/common/deferred.py:408: FutureWarning: `Value.case` is deprecated as of v10.0.0; use value.cases() or ibis.cases()\n",
+      "  return func(*args, **kwargs)\n"
+     ]
+    },
+    {
+     "ename": "NameError",
+     "evalue": "name 'non_conserved' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
+      "\u001b[31mNameError\u001b[39m                                 Traceback (most recent call last)",
+      "\u001b[36mFile \u001b[39m\u001b[32m<timed exec>:50\u001b[39m\n",
+      "\u001b[31mNameError\u001b[39m: name 'non_conserved' is not defined"
+     ]
+    }
+   ],
    "source": [
+    "%%time \n",
     "counties = con.read_parquet('../CA_counties.parquet')\n",
     "# ca = con.read_parquet(ca_temp_parquet)\n",
     "\n",
@@ -570,7 +660,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.8"
+   "version": "3.12.9"
   }
  },
  "nbformat": 4,

preprocess/utils.py

@@ -0,0 +1,198 @@
+from minio.error import S3Error
+
+import zipfile
+import os
+import subprocess
+
+import geopandas as gpd
+import ibis
+from ibis import _
+
+import rasterio
+from rasterio.features import shapes
+from shapely.geometry import shape
+import numpy as np
+
+
+def info(folder, file, bucket = "public-ca30x30", base_folder = 'CBN-data/'):
+    """
+    Extract minio path to upload/download data 
+    """
+    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):
+    """
+    Downloading file from minio
+    """
+    if not file_name: 
+        file_name = file
+    bucket, path = info(folder, file)
+    s3.fget_object(bucket, path ,file_name) 
+    return
+
+def upload(s3, folder, file):
+    """
+    Uploading file from minio
+    """
+    bucket, path = info(folder, file)
+    s3.fput_object(bucket, path ,file) 
+    return
+
+def unzip(folder, file):
+    """
+    Unzipping zip files 
+    """
+    download(s3, folder, file)
+    with zipfile.ZipFile(file, 'r') as zip_ref:
+        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"):
+    """
+    Driver function to process vectors 
+    """
+    if gdf is None:
+        gdf = gpd.read_file(file)
+    if gdf.crs != crs:
+        gdf = gdf.to_crs(crs)
+    if gdf.geometry.name != 'geom':
+        gdf = gdf.rename_geometry('geom')
+    if file_name:
+        file = file_name
+    # upload_parquet(folder, file, gdf)
+    name, ext = os.path.splitext(file)
+    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  
+
+
+def reproject_raster(input_file, crs="EPSG:3310"):
+    """
+    Reproject rasters 
+    """
+    suffix = '_processed'
+    name, ext = os.path.splitext(input_file)
+    output_file = f"{name}{suffix}{ext}"
+    command = [
+        "gdalwarp",
+        "-t_srs", crs,
+        input_file,
+        output_file 
+        ]
+    try:
+        subprocess.run(command, check=True)
+        print(f"Reprojection successful!")
+    except subprocess.CalledProcessError as e:
+        print(f"Error occurred during reprojection: {e}")
+    return output_file 
+
+def make_cog(input_file, crs="EPSG:4326"):
+    """
+    Converting TIF to COGs
+    """
+    suffix = '_COG'
+    name, ext = os.path.splitext(input_file)
+    output_file = f"{name}{suffix}{ext}"
+    command = [
+        "gdalwarp",
+        "-t_srs", crs,
+        "-of", "COG",
+        input_file,
+        output_file 
+        ]
+    try:
+        subprocess.run(command, check=True)
+        print(f"Successful!")
+    except subprocess.CalledProcessError as e:
+        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
+    
+        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)
+    print(gdf)
+    return output_file
+
+def filter_raster(s3, folder, file, percentile):
+    """
+    Helper function to filter rasteres 
+    """
+    with rasterio.open(file) as src:
+        data = src.read(1)  # Read the first band
+        profile = src.profile
+    # mask no data values
+    masked_data = np.ma.masked_equal(data, src.nodata)
+
+    # compute percentile/threshold 
+    p = np.percentile(masked_data.compressed(),percentile)
+    filtered = np.where(data >= p, data, src.nodata)
+    name, ext = os.path.splitext(file)
+    new_file = f"{name}{'_'}{percentile}{'percentile'}{ext}"
+
+    profile.update(dtype=rasterio.float64)
+    with rasterio.open(new_file, "w", **profile) as dst:
+        dst.write(filtered, 1)
+    process_raster(s3, folder, file)
+    return
+
+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):
+    """
+    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')
+    to_pmtiles(name+'.geojson', name+'.pmtiles', options = ['--extend-zooms-if-still-dropping'])
+    upload(s3, folder, name+'.pmtiles')
+    return
+