{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":64733,"databundleVersionId":7762802,"sourceType":"competition"},{"sourceId":7752763,"sourceType":"datasetVersion","datasetId":4533016}],"dockerImageVersionId":30664,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"Published on March 03, 2024. By Marília Prata, mpwolke","metadata":{}},{"cell_type":"markdown","source":"### In this notebook, various examples will be presented to preprocess and manipulate geographic data to create a species distribution model for predicting the probability of distribution.\n\n* Extracting values of desired locations from raster data\n* Filling missing values using spatial features\n* Merging data considering spatio-temporal information\n* Building a prediction model\n* Predicting and plotting distribution probability maps\n* Drawing variable importance maps\n* However, the provided image data in the competition will not be used. It seems that developing a multimodal model by incorporating image information in the future could lead to a high-performance model. Please note that the model presented here is a baseline model.","metadata":{}},{"cell_type":"code","source":"from glob import glob\nimport os\n# from netCDF4 import Dataset\nimport pandas as pd\nimport geopandas as gpd #\nfrom shapely.geometry import Polygon, LineString, Point\nimport rasterio\nfrom rasterio.plot import show\n# from rasterio.transform import from_origin\nfrom rasterstats import zonal_stats\nimport matplotlib.pyplot as plt\nimport cv2\nimport numpy as np\nimport geopandas as gpd\nimport tqdm\n\nroot_path = '/kaggle/input'","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:20.625511Z","iopub.execute_input":"2024-03-12T09:25:20.625807Z","iopub.status.idle":"2024-03-12T09:25:22.082395Z","shell.execute_reply.started":"2024-03-12T09:25:20.625785Z","shell.execute_reply":"2024-03-12T09:25:22.081407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Read rasterfile and Visualize with point vector","metadata":{}},{"cell_type":"code","source":"# Read rasterfile\nbio_raster_path_list = glob(root_path + '/climateclef/Climate/BioClimatic_Average_1981-2010/*.tif')\nclimatic_raster_path_list = glob(root_path + '/climateclef/Climate/Climatic_Monthly_2000-2019/*.tif')","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:22.083739Z","iopub.execute_input":"2024-03-12T09:25:22.084093Z","iopub.status.idle":"2024-03-12T09:25:22.366755Z","shell.execute_reply.started":"2024-03-12T09:25:22.084073Z","shell.execute_reply":"2024-03-12T09:25:22.366168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot raster file\nfor i,raster_path in enumerate(climatic_raster_path_list[:2]):\n    with rasterio.open(raster_path) as src:\n        plt.figure(figsize = [8,8])\n        plt.title(raster_path.split(\"/\")[-1])\n        show(src)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:22.367606Z","iopub.execute_input":"2024-03-12T09:25:22.367968Z","iopub.status.idle":"2024-03-12T09:25:28.656599Z","shell.execute_reply.started":"2024-03-12T09:25:22.367931Z","shell.execute_reply":"2024-03-12T09:25:28.655659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot raster file\nfor i,raster_path in enumerate(bio_raster_path_list[:2]):\n    with rasterio.open(raster_path) as src:\n        plt.figure(figsize = [8,8])\n        plt.title(raster_path.split(\"/\")[-1])\n        show(src)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:28.658975Z","iopub.execute_input":"2024-03-12T09:25:28.659314Z","iopub.status.idle":"2024-03-12T09:25:33.564745Z","shell.execute_reply.started":"2024-03-12T09:25:28.659286Z","shell.execute_reply":"2024-03-12T09:25:33.563861Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize every monitoring site in meta data\n# Example for 1000 samples\n# Read meta data\ntrain_meta = pd.read_csv(root_path + \"/geolifeclef-2024/GLC24_PA_metadata_train.csv\")\ntrain_meta.head()\n\nsub_train_meta = train_meta.drop_duplicates('surveyId').sample(n=1000, random_state=42)\nsub_train_meta.index = range(len(sub_train_meta))\n\n# make Point vector \npoint_list = []\nfor i in tqdm.tqdm(range(len(sub_train_meta))):\n    x,y = sub_train_meta.loc[i, ['lon', 'lat']]\n    poind_i = Point(x,y)\n    point_list.append(poind_i)\n\nsub_train_meta.loc[:,'geometry'] = point_list\n\n# Read meta data\ntest_meta = pd.read_csv(root_path + \"/geolifeclef-2024/GLC24_PA_metadata_test.csv\")\ntest_meta.head()\n\nsub_test_meta = test_meta.drop_duplicates('surveyId').sample(n=1000, random_state=42)\nsub_test_meta.index = range(len(sub_test_meta))\n\n# make Point vector \npoint_list = []\nfor i in tqdm.tqdm(range(len(sub_test_meta))):\n    x,y = sub_test_meta.loc[i, ['lon', 'lat']]\n    poind_i = Point(x,y)\n    point_list.append(poind_i)\n\nsub_test_meta.loc[:,'geometry'] = point_list","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:33.565825Z","iopub.execute_input":"2024-03-12T09:25:33.566025Z","iopub.status.idle":"2024-03-12T09:25:35.976841Z","shell.execute_reply.started":"2024-03-12T09:25:33.566005Z","shell.execute_reply":"2024-03-12T09:25:35.976174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sub_meta_gdf = gpd.GeoDataFrame(sub_train_meta, geometry = 'geometry')\ntrain_sub_meta_gdf.crs = {'init':'epsg:4326'} \ntrain_sub_meta_gdf.head()\n\ntest_sub_meta_gdf = gpd.GeoDataFrame(sub_test_meta, geometry = 'geometry')\ntest_sub_meta_gdf.crs = {'init':'epsg:4326'} \ntest_sub_meta_gdf.head()\n\n# visualize each monitoring site\nm = train_sub_meta_gdf.drop_duplicates(['lon', 'lat']).explore(color = 'green')\ntest_sub_meta_gdf.drop_duplicates(['lon', 'lat']).explore(m=m, color = 'red')","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:35.977828Z","iopub.execute_input":"2024-03-12T09:25:35.978046Z","iopub.status.idle":"2024-03-12T09:25:40.676371Z","shell.execute_reply.started":"2024-03-12T09:25:35.978025Z","shell.execute_reply":"2024-03-12T09:25:40.675477Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# visualize with raster file\nfig, ax = plt.subplots(figsize=(8, 8))\n\nraster_path = bio_raster_path_list[0]\nwith rasterio.open(raster_path) as src:\n    show(src, ax=ax)\ntrain_sub_meta_gdf.drop_duplicates(['lon', 'lat']).plot(ax=ax,color = 'green')\ntest_sub_meta_gdf.drop_duplicates(['lon', 'lat']).plot(ax=ax,color = 'red')","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:40.677356Z","iopub.execute_input":"2024-03-12T09:25:40.677709Z","iopub.status.idle":"2024-03-12T09:25:45.326959Z","shell.execute_reply.started":"2024-03-12T09:25:40.677689Z","shell.execute_reply":"2024-03-12T09:25:45.326058Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:45.327909Z","iopub.execute_input":"2024-03-12T09:25:45.328214Z","iopub.status.idle":"2024-03-12T09:25:45.346275Z","shell.execute_reply.started":"2024-03-12T09:25:45.328186Z","shell.execute_reply":"2024-03-12T09:25:45.345349Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### merge raster data with meta data\nThis competition provides a CSV file by completing the tasks below. However, I have written example code below for practicing data preprocessing on your own","metadata":{}},{"cell_type":"code","source":"%%time\n# two raster file example\ntrain_geo_info_gdf = train_sub_meta_gdf.drop_duplicates(['lon', 'lat'])\ntrain_geo_info_gdf.index = range(len(train_geo_info_gdf))\n\ntest_geo_info_gdf = test_sub_meta_gdf.drop_duplicates(['lon', 'lat'])\ntest_geo_info_gdf.index = range(len(test_geo_info_gdf))\n\nfor raster_path in tqdm.tqdm(bio_raster_path_list[:2]):\n    bio_raster = rasterio.open(raster_path)\n    \n    bio_values = zonal_stats(vectors=train_geo_info_gdf.geometry,raster=bio_raster.read(1), affine=bio_raster.transform,stats=['mean'])\n    bio_values_list = [value['mean'] for value in bio_values]\n    feature_name = raster_path.split(\"/\")[-1].split(\".\")[0]\n    train_geo_info_gdf.loc[:,feature_name] = bio_values_list\n    \n    bio_values = zonal_stats(vectors=test_geo_info_gdf.geometry,raster=bio_raster.read(1), affine=bio_raster.transform,stats=['mean'])\n    bio_values_list = [value['mean'] for value in bio_values]\n    feature_name = raster_path.split(\"/\")[-1].split(\".\")[0]\n    test_geo_info_gdf.loc[:,feature_name] = bio_values_list","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:45.347606Z","iopub.execute_input":"2024-03-12T09:25:45.347926Z","iopub.status.idle":"2024-03-12T09:25:50.213832Z","shell.execute_reply.started":"2024-03-12T09:25:45.347887Z","shell.execute_reply":"2024-03-12T09:25:50.213223Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_geo_info_gdf.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:50.216480Z","iopub.execute_input":"2024-03-12T09:25:50.217137Z","iopub.status.idle":"2024-03-12T09:25:50.232852Z","shell.execute_reply.started":"2024-03-12T09:25:50.217090Z","shell.execute_reply":"2024-03-12T09:25:50.231991Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_geo_info_gdf.head()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:50.233983Z","iopub.execute_input":"2024-03-12T09:25:50.234257Z","iopub.status.idle":"2024-03-12T09:25:50.257832Z","shell.execute_reply.started":"2024-03-12T09:25:50.234233Z","shell.execute_reply":"2024-03-12T09:25:50.256943Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Dataset merge\nUtilize the CSV file provided in this competition\n\n\n* landsat timeseries data has missing value\n* If Landsat timeseries data is collected after the sampling date, it is considered as data not available for the originally intended prediction timeframe. If not treated as missing values, it may lead to data leakage issues.\n* We need to consider this timeframe issue when merging datasets.","metadata":{}},{"cell_type":"markdown","source":"* make train dataset","metadata":{}},{"cell_type":"code","source":"train_meta = pd.read_csv(root_path + \"/geolifeclef-2024/GLC24_PA_metadata_train.csv\")","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:50.258975Z","iopub.execute_input":"2024-03-12T09:25:50.259233Z","iopub.status.idle":"2024-03-12T09:25:51.083278Z","shell.execute_reply.started":"2024-03-12T09:25:50.259212Z","shell.execute_reply":"2024-03-12T09:25:51.082417Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# read average landsat dataset\n# Simply merge; this set has no missing values.\ntrain_bioclimatic = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Climate/Average 1981-2010/GLC24-PA-train-bioclimatic.csv\")\ntrain_bioclimatic_monthly = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Climate/Monthly/GLC24-PA-train-bioclimatic_monthly.csv\")\ntrain_elevation = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Elevation/GLC24-PA-train-elevation.csv\")\ntrain_human_footprint = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Human Footprint/GLC24-PA-train-human_footprint.csv\")\ntrain_land_cover = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/LandCover/GLC24-PA-train-landcover.csv\")\ntrain_soilgrid = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/SoilGrids/GLC24-PA-train-soilgrids.csv\")\n\n# first merge bioclimatic averate dataset\nmerge_key = 'surveyId'\n\n# merge tables\ntrain_merged = pd.merge(train_bioclimatic, train_bioclimatic_monthly, on=merge_key)\ntrain_merged = pd.merge(train_merged, train_elevation, on=merge_key)\ntrain_merged = pd.merge(train_merged, train_human_footprint, on=merge_key)\ntrain_merged = pd.merge(train_merged, train_land_cover, on=merge_key)\ntrain_merged = pd.merge(train_merged, train_soilgrid, on=merge_key)\n\ntrain_merged","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:25:51.084360Z","iopub.execute_input":"2024-03-12T09:25:51.084608Z","iopub.status.idle":"2024-03-12T09:26:00.360684Z","shell.execute_reply.started":"2024-03-12T09:25:51.084586Z","shell.execute_reply":"2024-03-12T09:26:00.359896Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# landsat time series missing value check\nlandsat_example = pd.read_csv(root_path + f\"/geolifeclef-2024/PA-train-landsat_time_series/GLC24-PA-train-landsat_time_series-blue.csv\")\nlandsat_example.head() ## missing value!","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:26:00.361808Z","iopub.execute_input":"2024-03-12T09:26:00.362116Z","iopub.status.idle":"2024-03-12T09:26:01.036091Z","shell.execute_reply.started":"2024-03-12T09:26:00.362073Z","shell.execute_reply":"2024-03-12T09:26:01.035297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We merge climatic timeseries data within a 10-year time window.(4 q * 10 years = 40)\n# For example, if the survey year was 2018, we merge monthly data from the years 2008 to 2017.\n\ndef shift_and_slice(array):\n    for i in range(len(array)):\n        row = array[i]\n        nan_indices = np.isnan(row)\n\n        rightmost_nan_index = np.argmax(nan_indices)\n\n        if np.any(nan_indices):\n            shifted_row = np.roll(row, len(row) - rightmost_nan_index)\n            array[i] = shifted_row\n    \n    return array[:,-(4 * 10):]\n\nfor monthly_landsat_data in tqdm.tqdm(['blue', 'green', 'red', 'nir', 'swir1', 'swir2']):\n    landsat_timeseries = pd.read_csv(root_path + f\"/geolifeclef-2024/PA-train-landsat_time_series/GLC24-PA-train-landsat_time_series-{monthly_landsat_data}.csv\")\n    origin_arr = landsat_timeseries.iloc[:,1:].values\n    new_arr = shift_and_slice(origin_arr)\n    \n    window_dataset =  landsat_timeseries.iloc[:,0:1]\n    window_dataset.loc[:,[f'{monthly_landsat_data}_t-{i+1}' for i in range(40)]] = new_arr\n    \n    train_merged = pd.merge(train_merged, window_dataset, on=merge_key)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:26:01.037158Z","iopub.execute_input":"2024-03-12T09:26:01.037385Z","iopub.status.idle":"2024-03-12T09:26:15.306652Z","shell.execute_reply.started":"2024-03-12T09:26:01.037365Z","shell.execute_reply":"2024-03-12T09:26:15.305416Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_merged","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:26:15.307769Z","iopub.execute_input":"2024-03-12T09:26:15.308082Z","iopub.status.idle":"2024-03-12T09:26:15.338855Z","shell.execute_reply.started":"2024-03-12T09:26:15.308056Z","shell.execute_reply":"2024-03-12T09:26:15.337979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\n# species id range from 1 to 11254\nprint(np.max(train_meta.drop_duplicates('speciesId').loc[:,['speciesId']].values), np.min(train_meta.drop_duplicates('speciesId').loc[:,['speciesId']].values))\nspecies_target_columns = [str(sp_id) for sp_id in range(1,11255)]\ninput_columns = train_merged.columns.tolist()[1:]\n\n# reshape metadata \nsurveyId_list = train_meta.drop_duplicates('surveyId').loc[:,['surveyId']].values\ntarget_arr = np.zeros([len(surveyId_list),11254])\nfor i,survey_id in tqdm.tqdm(enumerate(surveyId_list)):\n    presense_species_ids = train_meta.loc[train_meta.surveyId == survey_id[0],'speciesId'].values\n    target_arr[i,[int(ids)-1 for ids in presense_species_ids]] = 1\n\ntrain_target_df = pd.DataFrame(target_arr)\ntrain_target_df.columns = species_target_columns\ntrain_target_df.loc[:,'surveyId'] = surveyId_list\ntrain_target_df","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:26:15.339877Z","iopub.execute_input":"2024-03-12T09:26:15.340237Z","iopub.status.idle":"2024-03-12T09:28:34.908301Z","shell.execute_reply.started":"2024-03-12T09:26:15.340208Z","shell.execute_reply":"2024-03-12T09:28:34.907435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* spatial interpolation for missing values","metadata":{}},{"cell_type":"code","source":"# check missing values\nmissing_value_col_list = []\nrows_with_missing_values = {}\nfor column in train_merged.columns:\n    missing_rows = train_merged.index[train_merged[column].isnull()].tolist()\n    if len(missing_rows) > 0:\n        rows_with_missing_values[column] = len(missing_rows)\n        missing_value_col_list.append(column)\n\nprint(rows_with_missing_values)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:34.909521Z","iopub.execute_input":"2024-03-12T09:28:34.909804Z","iopub.status.idle":"2024-03-12T09:28:35.079482Z","shell.execute_reply.started":"2024-03-12T09:28:34.909777Z","shell.execute_reply":"2024-03-12T09:28:35.078665Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check missing value region\ntarget_colname = 'swir2_t-39'\nsub_df = pd.merge(train_merged.loc[:, ['surveyId', target_colname]], train_meta.loc[:, ['lon', 'lat', 'surveyId']].drop_duplicates(\"surveyId\"), on='surveyId')\nna_ind = sub_df.loc[:,target_colname].isna().values\n\nsub_df_train = sub_df.loc[~na_ind,:]\nsub_df_target = sub_df.loc[na_ind,:]\nx,y,z = sub_df_train['lon'].values, sub_df_train['lat'].values, sub_df_train[target_colname].values\nx_target,y_target,z_target = sub_df_target['lon'].values, sub_df_target['lat'].values, np.zeros(sub_df_target[target_colname].shape)\n\nplt.figure(figsize=(12, 6))\nplt.scatter(x, y, marker='o', color='green', label='Non Missing Data')\nplt.scatter(x_target, y_target, marker='o', color='red', label='Missing Data', s = 5)\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:35.080564Z","iopub.execute_input":"2024-03-12T09:28:35.081064Z","iopub.status.idle":"2024-03-12T09:28:35.796654Z","shell.execute_reply.started":"2024-03-12T09:28:35.081043Z","shell.execute_reply":"2024-03-12T09:28:35.795971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check missing value region\ntarget_colname = 'Elevation'\nsub_df = pd.merge(train_merged.loc[:, ['surveyId', target_colname]], train_meta.loc[:, ['lon', 'lat', 'surveyId']].drop_duplicates(\"surveyId\"), on='surveyId')\nna_ind = sub_df.loc[:,target_colname].isna().values\n\nsub_df_train = sub_df.loc[~na_ind,:]\nsub_df_target = sub_df.loc[na_ind,:]\nx,y,z = sub_df_train['lon'].values, sub_df_train['lat'].values, sub_df_train[target_colname].values\nx_target,y_target,z_target = sub_df_target['lon'].values, sub_df_target['lat'].values, np.zeros(sub_df_target[target_colname].shape)\n\nplt.figure(figsize=(12, 6))\nplt.scatter(x, y, marker='o', color='green', label='Non Missing Data')\nplt.scatter(x_target, y_target, marker='o', color='red', label='Missing Data')\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:35.797813Z","iopub.execute_input":"2024-03-12T09:28:35.798270Z","iopub.status.idle":"2024-03-12T09:28:36.504837Z","shell.execute_reply.started":"2024-03-12T09:28:35.798242Z","shell.execute_reply":"2024-03-12T09:28:36.503759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check missing value region\ntarget_colname = 'HumanFootprint-Roads'\nsub_df = pd.merge(train_merged.loc[:, ['surveyId', target_colname]], train_meta.loc[:, ['lon', 'lat', 'surveyId']].drop_duplicates(\"surveyId\"), on='surveyId')\nna_ind = sub_df.loc[:,target_colname].isna().values\n\nsub_df_train = sub_df.loc[~na_ind,:]\nsub_df_target = sub_df.loc[na_ind,:]\nx,y,z = sub_df_train['lon'].values, sub_df_train['lat'].values, sub_df_train[target_colname].values\nx_target,y_target,z_target = sub_df_target['lon'].values, sub_df_target['lat'].values, np.zeros(sub_df_target[target_colname].shape)\n\nplt.figure(figsize=(12, 6))\nplt.scatter(x, y, marker='o', color='green', label='Non Missing Data')\nplt.scatter(x_target, y_target, marker='o', color='red', label='Missing Data', s = 5)\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:36.506434Z","iopub.execute_input":"2024-03-12T09:28:36.506967Z","iopub.status.idle":"2024-03-12T09:28:37.219929Z","shell.execute_reply.started":"2024-03-12T09:28:36.506924Z","shell.execute_reply":"2024-03-12T09:28:37.219040Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check missing value region\ntarget_colname = 'Soilgrid-bdod'\nsub_df = pd.merge(train_merged.loc[:, ['surveyId', target_colname]], train_meta.loc[:, ['lon', 'lat', 'surveyId']].drop_duplicates(\"surveyId\"), on='surveyId')\nna_ind = sub_df.loc[:,target_colname].isna().values\n\nsub_df_train = sub_df.loc[~na_ind,:]\nsub_df_target = sub_df.loc[na_ind,:]\nx,y,z = sub_df_train['lon'].values, sub_df_train['lat'].values, sub_df_train[target_colname].values\nx_target,y_target,z_target = sub_df_target['lon'].values, sub_df_target['lat'].values, np.zeros(sub_df_target[target_colname].shape)\n\nplt.figure(figsize=(12, 6))\nplt.scatter(x, y, marker='o', color='green', label='Non Missing Data')\nplt.scatter(x_target, y_target, marker='o', color='red', label='Missing Data', s = 5)\n\nplt.legend()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:37.220831Z","iopub.execute_input":"2024-03-12T09:28:37.221050Z","iopub.status.idle":"2024-03-12T09:28:37.954879Z","shell.execute_reply.started":"2024-03-12T09:28:37.221030Z","shell.execute_reply":"2024-03-12T09:28:37.954317Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Based on observations from the plot, it appears that there are densely populated measurements near the missing regions. \n* To handle interpolation, a straightforward approach involves using the values from the nearest neighbors. If a more reliable method is desired, considering techniques such as kriging is also worth exploring.","metadata":{}},{"cell_type":"code","source":"%%time\n# missing_value_col_list\nfor target_colname in tqdm.tqdm(missing_value_col_list):\n\n    sub_df = pd.merge(train_merged.loc[:, ['surveyId', target_colname]], train_meta.loc[:, ['lon', 'lat', 'surveyId']].drop_duplicates(\"surveyId\"), on='surveyId')\n    na_ind_bool = sub_df.loc[:,target_colname].isna().values\n\n    na_ind = sub_df.index[na_ind_bool]\n    no_na_ind = sub_df.index[~na_ind_bool]\n\n\n\n    given_points = sub_df.iloc[no_na_ind,[2,3]].values\n\n    for na_ind_i in na_ind:\n\n        target_point = sub_df.iloc[na_ind_i,[2,3]].values\n        distances = np.linalg.norm(given_points - target_point, axis=1)\n        closest_index = no_na_ind[np.argmin(distances)]\n        closest_value = sub_df.iloc[closest_index,1]\n        train_merged.loc[na_ind_i, target_colname] = closest_value\n    ","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:28:37.955934Z","iopub.execute_input":"2024-03-12T09:28:37.956322Z","iopub.status.idle":"2024-03-12T09:29:53.484064Z","shell.execute_reply.started":"2024-03-12T09:28:37.956301Z","shell.execute_reply":"2024-03-12T09:29:53.483145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check missing values\nrows_with_missing_values = {}\nfor column in train_merged.columns:\n    missing_rows = train_merged.index[train_merged[column].isnull()].tolist()\n    if len(missing_rows) > 0:\n        rows_with_missing_values[column] = len(missing_rows)\nprint(rows_with_missing_values)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:29:53.484998Z","iopub.execute_input":"2024-03-12T09:29:53.485221Z","iopub.status.idle":"2024-03-12T09:29:53.652031Z","shell.execute_reply.started":"2024-03-12T09:29:53.485200Z","shell.execute_reply":"2024-03-12T09:29:53.651156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check -inf values\ncolumns_with_negative_inf = train_merged.columns[train_merged.apply(lambda col: any(col == float('-inf')))]\nif len(columns_with_negative_inf) > 0:\n    print(\"This column contains a value of -inf\")\n    for column_name in columns_with_negative_inf:\n        print(column_name)\nelse:\n    print(\"no inf value\")\n","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:29:53.652838Z","iopub.execute_input":"2024-03-12T09:29:53.653066Z","iopub.status.idle":"2024-03-12T09:30:02.304273Z","shell.execute_reply.started":"2024-03-12T09:29:53.653047Z","shell.execute_reply":"2024-03-12T09:30:02.303493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Replace -inf with zero.\ntrain_merged.replace({float('-inf'): 0}, inplace=True)\n\n# check -inf values\ncolumns_with_negative_inf = train_merged.columns[train_merged.apply(lambda col: any(col == float('-inf')))]\nif len(columns_with_negative_inf) > 0:\n    print(\"This column contains a value of -inf\")\n    for column_name in columns_with_negative_inf:\n        print(column_name)\nelse:\n    print(\"no inf value\")","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:30:02.305422Z","iopub.execute_input":"2024-03-12T09:30:02.306328Z","iopub.status.idle":"2024-03-12T09:30:11.022143Z","shell.execute_reply.started":"2024-03-12T09:30:02.306300Z","shell.execute_reply":"2024-03-12T09:30:11.021331Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* make test dataset","metadata":{}},{"cell_type":"code","source":"%%time\n# read average landsat dataset\n# Simply merge; this set has no missing values.\ntest_bioclimatic = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Climate/Average 1981-2010/GLC24-PA-test-bioclimatic.csv\")\ntest_bioclimatic_monthly = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Climate/Monthly/GLC24-PA-test-bioclimatic_monthly.csv\")\ntest_elevation = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Elevation/GLC24-PA-test-elevation.csv\")\ntest_human_footprint = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/Human Footprint/GLC24-PA-test-human_footprint.csv\")\ntest_land_cover = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/LandCover/GLC24-PA-test-landcover.csv\")\ntest_soilgrid = pd.read_csv(root_path + \"/geolifeclef-2024/EnvironmentalRasters/EnvironmentalRasters/SoilGrids/GLC24-PA-test-soilgrids.csv\")\n\n# first merge bioclimatic averate dataset\nmerge_key = 'surveyId'\n\n# merge tables\ntest_merged = pd.merge(test_bioclimatic, test_bioclimatic_monthly, on=merge_key)\ntest_merged = pd.merge(test_merged, test_elevation, on=merge_key)\ntest_merged = pd.merge(test_merged, test_human_footprint, on=merge_key)\ntest_merged = pd.merge(test_merged, test_land_cover, on=merge_key)\ntest_merged = pd.merge(test_merged, test_soilgrid, on=merge_key)\n\ntest_merged","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:30:11.023170Z","iopub.execute_input":"2024-03-12T09:30:11.023439Z","iopub.status.idle":"2024-03-12T09:30:11.628335Z","shell.execute_reply.started":"2024-03-12T09:30:11.023414Z","shell.execute_reply":"2024-03-12T09:30:11.627498Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# We merge climatic timeseries data within a 10-year time window.(4 q * 10 years = 40)\n# For example, if the survey year was 2018, we merge monthly data from the years 2008 to 2017.\nfor monthly_landsat_data in tqdm.tqdm(['blue', 'green', 'red', 'nir', 'swir1', 'swir2']):\n    landsat_timeseries = pd.read_csv(root_path + f\"/geolifeclef-2024/PA-test-landsat_time_series/GLC24-PA-test-landsat_time_series-{monthly_landsat_data}.csv\")\n    origin_arr = landsat_timeseries.iloc[:,1:].values\n    new_arr = shift_and_slice(origin_arr)\n    \n    window_dataset =  landsat_timeseries.iloc[:,0:1]\n    window_dataset.loc[:,[f'{monthly_landsat_data}_t-{i+1}' for i in range(40)]] = new_arr\n    \n    test_merged = pd.merge(test_merged, window_dataset, on=merge_key)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:30:11.632253Z","iopub.execute_input":"2024-03-12T09:30:11.632497Z","iopub.status.idle":"2024-03-12T09:30:12.524600Z","shell.execute_reply.started":"2024-03-12T09:30:11.632476Z","shell.execute_reply":"2024-03-12T09:30:12.523834Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_merged","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:30:12.525533Z","iopub.execute_input":"2024-03-12T09:30:12.525753Z","iopub.status.idle":"2024-03-12T09:30:12.551133Z","shell.execute_reply.started":"2024-03-12T09:30:12.525734Z","shell.execute_reply":"2024-03-12T09:30:12.550178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* train and prediction example for speciesId 540","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nfrom xgboost import XGBClassifier\nimport xgboost as xgb","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:30:12.552264Z","iopub.execute_input":"2024-03-12T09:30:12.552650Z","iopub.status.idle":"2024-03-12T09:30:12.706478Z","shell.execute_reply.started":"2024-03-12T09:30:12.552626Z","shell.execute_reply":"2024-03-12T09:30:12.705667Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train\ntartget_speciesId = '540'\n\ntrain_df_sub = pd.merge(train_target_df.loc[:,['surveyId', tartget_speciesId]], train_merged, on='surveyId')\ntrain_df_sub.head()\n\nX, y = train_df_sub.iloc[:,2:].values, train_df_sub.iloc[:,1].values\nX_train, X_valid, y_train, y_valid = train_test_split(X, y, test_size=0.1, random_state=7)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:39:59.405417Z","iopub.execute_input":"2024-03-12T09:39:59.405751Z","iopub.status.idle":"2024-03-12T09:40:00.505228Z","shell.execute_reply.started":"2024-03-12T09:39:59.405726Z","shell.execute_reply":"2024-03-12T09:40:00.504353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.hist(y)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:40:00.506628Z","iopub.execute_input":"2024-03-12T09:40:00.506907Z","iopub.status.idle":"2024-03-12T09:40:00.692232Z","shell.execute_reply.started":"2024-03-12T09:40:00.506886Z","shell.execute_reply":"2024-03-12T09:40:00.691443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\nmodel = XGBClassifier(objective='binary:logistic',eval_metric=['error','logloss'], verbose = 1)\nmodel.fit(X_train, y_train, eval_set=[(X_valid, y_valid)], verbose=True,)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:40:12.786985Z","iopub.execute_input":"2024-03-12T09:40:12.787345Z","iopub.status.idle":"2024-03-12T09:41:06.047903Z","shell.execute_reply.started":"2024-03-12T09:40:12.787319Z","shell.execute_reply":"2024-03-12T09:41:06.047138Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_merged","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:41:06.049415Z","iopub.execute_input":"2024-03-12T09:41:06.049740Z","iopub.status.idle":"2024-03-12T09:41:06.075640Z","shell.execute_reply.started":"2024-03-12T09:41:06.049711Z","shell.execute_reply":"2024-03-12T09:41:06.074818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* prediction and mapping species distributions","metadata":{}},{"cell_type":"code","source":"test_meta = pd.read_csv(root_path + \"/geolifeclef-2024/GLC24_PA_metadata_test.csv\")\nX_test = test_merged.iloc[:,1:].values\ny_pred_proba = model.predict_proba(X_test)[:,1]\ntest_meta.loc[:,tartget_speciesId + \"_prob\"] = y_pred_proba","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:08.032400Z","iopub.execute_input":"2024-03-12T09:53:08.032724Z","iopub.status.idle":"2024-03-12T09:53:08.081436Z","shell.execute_reply.started":"2024-03-12T09:53:08.032701Z","shell.execute_reply":"2024-03-12T09:53:08.080570Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:08.455935Z","iopub.execute_input":"2024-03-12T09:53:08.456805Z","iopub.status.idle":"2024-03-12T09:53:08.471186Z","shell.execute_reply.started":"2024-03-12T09:53:08.456779Z","shell.execute_reply":"2024-03-12T09:53:08.470381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# make Point vector \npoint_list = []\nfor i in tqdm.tqdm(range(len(test_meta))):\n    x,y = test_meta.loc[i, ['lon', 'lat']]\n    poind_i = Point(x,y)\n    point_list.append(poind_i)\n\ntest_meta.loc[:,'geometry'] = point_list","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:08.985880Z","iopub.execute_input":"2024-03-12T09:53:08.986234Z","iopub.status.idle":"2024-03-12T09:53:10.855381Z","shell.execute_reply.started":"2024-03-12T09:53:08.986208Z","shell.execute_reply":"2024-03-12T09:53:10.854352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta_gdf = gpd.GeoDataFrame(test_meta, geometry = 'geometry')\ntest_meta_gdf.crs = {'init':'epsg:4326'} ","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:10.857536Z","iopub.execute_input":"2024-03-12T09:53:10.857893Z","iopub.status.idle":"2024-03-12T09:53:10.903695Z","shell.execute_reply.started":"2024-03-12T09:53:10.857863Z","shell.execute_reply":"2024-03-12T09:53:10.902763Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta_gdf.explore(column = tartget_speciesId + \"_prob\", cmap = 'jet', legend = True)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:10.904782Z","iopub.execute_input":"2024-03-12T09:53:10.905001Z","iopub.status.idle":"2024-03-12T09:53:12.022369Z","shell.execute_reply.started":"2024-03-12T09:53:10.904981Z","shell.execute_reply":"2024-03-12T09:53:12.021521Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* feature importance","metadata":{}},{"cell_type":"code","source":"import shap\n\nexplainer = shap.TreeExplainer(model)\nshap_values = explainer(X_test)\nfeatures_col_list = train_merged.columns[1:]\n\nshap.summary_plot(shap_values, X_test, feature_names = features_col_list)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:14.702008Z","iopub.execute_input":"2024-03-12T09:53:14.702344Z","iopub.status.idle":"2024-03-12T09:53:18.297627Z","shell.execute_reply.started":"2024-03-12T09:53:14.702319Z","shell.execute_reply":"2024-03-12T09:53:18.296802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* feature importance map Bio-pr_04_2013","metadata":{}},{"cell_type":"code","source":"target_feature = 'Bio-pr_04_2013'\nfeature_shap_values = shap_values[:,features_col_list == target_feature].values","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:49.244132Z","iopub.execute_input":"2024-03-12T09:53:49.244435Z","iopub.status.idle":"2024-03-12T09:53:49.297032Z","shell.execute_reply.started":"2024-03-12T09:53:49.244413Z","shell.execute_reply":"2024-03-12T09:53:49.296310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_shap_values.shape","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:49.728397Z","iopub.execute_input":"2024-03-12T09:53:49.728715Z","iopub.status.idle":"2024-03-12T09:53:49.734939Z","shell.execute_reply.started":"2024-03-12T09:53:49.728693Z","shell.execute_reply":"2024-03-12T09:53:49.733847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta_gdf.loc[:,'feature_contribution'] = feature_shap_values","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:53:52.547382Z","iopub.execute_input":"2024-03-12T09:53:52.548013Z","iopub.status.idle":"2024-03-12T09:53:52.553357Z","shell.execute_reply.started":"2024-03-12T09:53:52.547981Z","shell.execute_reply":"2024-03-12T09:53:52.552448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_meta_gdf.explore(column = 'feature_contribution', cmap = 'jet', legend = True)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:54:08.361373Z","iopub.execute_input":"2024-03-12T09:54:08.361694Z","iopub.status.idle":"2024-03-12T09:54:09.531696Z","shell.execute_reply.started":"2024-03-12T09:54:08.361672Z","shell.execute_reply":"2024-03-12T09:54:09.530629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* feature importance map elevation","metadata":{}},{"cell_type":"code","source":"target_feature = 'Elevation'\nfeature_shap_values = shap_values[:,features_col_list == target_feature].values\ntest_meta_gdf.loc[:,'feature_contribution'] = feature_shap_values\ntest_meta_gdf.explore(column = 'feature_contribution', cmap = 'jet', legend = True)","metadata":{"execution":{"iopub.status.busy":"2024-03-12T09:55:17.447300Z","iopub.execute_input":"2024-03-12T09:55:17.447632Z","iopub.status.idle":"2024-03-12T09:55:19.168048Z","shell.execute_reply.started":"2024-03-12T09:55:17.447609Z","shell.execute_reply":"2024-03-12T09:55:19.166857Z"},"trusted":true},"execution_count":null,"outputs":[]}]}