<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 0. Standard setup for logging and plotting inside a notebook
Step1:  1. Choose a representative species for a case study
Step2:  2. Rasterize the species, to get a matrix of pixels
Step3:  2.1 Plot to get an idea
Step4:  3. Load the biogeographical regons raster layer
Step5:  3.1 Plot to get an idea
Step6:  3.2 Load the continents vector layer (for further clipping of pseudo-absence area), rasterize
Step7:  As agreed, we will merge Europe and Asia to be a bit closer to the biogeographical regions. We do this user-specific patching for now, until a better solution is found.
Step8:  Set all values >1 to 1. (we only care about presence/absence)
Step9:  Delete band 2 (Europe, previously merged with layer 0==Asia)
Step10:  This is how the band 0 looks like now.
Step11:  4. Sample pseudo-absence pixels, taking into account all the distinct biomes that fall in the species region.
Step12:  4.1 Plot the biomes taken into account for sampling pseudo-absences, to get an idea
Step13:  4.2 Plot the sampled pseudo-absences, to get an idea
Step14:  5. Construct a convenient dataframe for testing with different SDM models
For the Example 2 datasheet, all cells of a global raster map are needed, one pixel per row.
5.1 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a "base" (zeroes) raster map.
Step15:  5.2 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a presences pixel map
Step16:  5.3 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a pseudo_absences pixel map
Step17:  5.4 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a minimum temperature pixel map
Step18:  5.5 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a maximum temperature pixel map
Step19:  5.6 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a mean temperature pixel map
Step20:  http
Step21:  6. Repeat with other species
Step22:  Third species.... (largest IUCN area, also plenty of occurrences in GBIF)
Step23:  Hmm, why does it take South Afrika into account?
Step24:  This above adds two bands with ones and zeroes
Step25:  Add continents column
<ASSISTANT_TASK:>
Python Code:
import logging
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
root = logging.getLogger()
root.addHandler(logging.StreamHandler())
%matplotlib inline
Explanation: 0. Standard setup for logging and plotting inside a notebook
End of explanation
# download from Google Drive: https://drive.google.com/open?id=0B9cazFzBtPuCOFNiUHYwcVFVODQ
# Representative example with multiple polygons in the shapefile, and a lot of point-records (also outside rangemaps)
from iSDM.species import IUCNSpecies
salmo_trutta = IUCNSpecies(name_species='Salmo trutta')
salmo_trutta.load_shapefile("../data/fish/selection/salmo_trutta")
Explanation: 1. Choose a representative species for a case study
End of explanation
rasterized = salmo_trutta.rasterize(raster_file="./salmo_trutta_full.tif", pixel_size=0.5)
Explanation: 2. Rasterize the species, to get a matrix of pixels
End of explanation
plt.figure(figsize=(25,20))
plt.imshow(rasterized, cmap="hot", interpolation="none")
Explanation: 2.1 Plot to get an idea
End of explanation
from iSDM.environment import RasterEnvironmentalLayer
biomes_adf = RasterEnvironmentalLayer(file_path="../data/rebioms/w001001.adf", name_layer="Biomes")
biomes_adf.load_data()
Explanation: 3. Load the biogeographical regons raster layer
End of explanation
biomes_adf.plot()
Explanation: 3.1 Plot to get an idea
End of explanation
from iSDM.environment import ContinentsLayer
from iSDM.environment import Source
continents = ContinentsLayer(file_path="../data/continents/continent.shp", source=Source.ARCGIS)
continents.load_data()
fig, ax = plt.subplots(1,1, figsize=(30,20))
continents.data_full.plot(column="continent", colormap="hsv")
continents_rasters = continents.rasterize(raster_file="../data/continents/continents_raster.tif", pixel_size=0.5)
continents_rasters.shape # stacked raster with 8 bands, one for each continent.
Explanation: 3.2 Load the continents vector layer (for further clipping of pseudo-absence area), rasterize
End of explanation
continents_rasters[0] = continents_rasters[0] + continents_rasters[2] # combine Europe and Asia
continents_rasters[0].max() # where the continents touch, we have overlap! that's why max is not 1, but 2.
Explanation: As agreed, we will merge Europe and Asia to be a bit closer to the biogeographical regions. We do this user-specific patching for now, until a better solution is found.
End of explanation
continents_rasters[0][continents_rasters[0] > 1] = 1
Explanation: Set all values >1 to 1. (we only care about presence/absence)
End of explanation
continents_rasters = np.delete(continents_rasters, 2, 0)
Explanation: Delete band 2 (Europe, previously merged with layer 0==Asia)
End of explanation
plt.figure(figsize=(25,20))
plt.imshow(continents_rasters[0], cmap="hot", interpolation="none")
continents_rasters.shape # now total of 7 band rather than 8
Explanation: This is how the band 0 looks like now.
End of explanation
selected_layers, pseudo_absences = biomes_adf.sample_pseudo_absences(species_raster_data=rasterized,continents_raster_data=continents_rasters, number_of_pseudopoints=1000)
Explanation: 4. Sample pseudo-absence pixels, taking into account all the distinct biomes that fall in the species region.
End of explanation
plt.figure(figsize=(25,20))
plt.imshow(selected_layers, cmap="hot", interpolation="none")
Explanation: 4.1 Plot the biomes taken into account for sampling pseudo-absences, to get an idea
End of explanation
plt.figure(figsize=(25,20))
plt.imshow(pseudo_absences, cmap="hot", interpolation="none")
Explanation: 4.2 Plot the sampled pseudo-absences, to get an idea
End of explanation
all_coordinates = biomes_adf.pixel_to_world_coordinates(raster_data=np.zeros_like(rasterized), filter_no_data_value=False)
all_coordinates
base_dataframe = pd.DataFrame([all_coordinates[0], all_coordinates[1]]).T
base_dataframe.columns=['decimallatitude', 'decimallongitude']
base_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
base_dataframe.head()
base_dataframe.tail()
Explanation: 5. Construct a convenient dataframe for testing with different SDM models
For the Example 2 datasheet, all cells of a global raster map are needed, one pixel per row.
5.1 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a "base" (zeroes) raster map.
End of explanation
presence_coordinates = salmo_trutta.pixel_to_world_coordinates()
presence_coordinates
presences_dataframe = pd.DataFrame([presence_coordinates[0], presence_coordinates[1]]).T
presences_dataframe.columns=['decimallatitude', 'decimallongitude']
presences_dataframe[salmo_trutta.name_species] = 1 # fill presences with 1's
presences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
presences_dataframe.head()
presences_dataframe.tail()
Explanation: 5.2 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a presences pixel map
End of explanation
pseudo_absence_coordinates = biomes_adf.pixel_to_world_coordinates(raster_data=pseudo_absences)
pseudo_absences_dataframe = pd.DataFrame([pseudo_absence_coordinates[0], pseudo_absence_coordinates[1]]).T
pseudo_absences_dataframe.columns=['decimallatitude', 'decimallongitude']
pseudo_absences_dataframe[salmo_trutta.name_species] = 0 # fill pseudo-absences with 0
pseudo_absences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
pseudo_absences_dataframe.head()
pseudo_absences_dataframe.tail()
Explanation: 5.3 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a pseudo_absences pixel map
End of explanation
from iSDM.environment import ClimateLayer
water_min_layer =  ClimateLayer(file_path="../data/watertemp/min_wt_2000.tif") 
water_min_reader = water_min_layer.load_data()
# HERE: should we ignore cells with no-data values for temperature? They are set to a really big negative number
# for now we keep them, otherwise could be NaN
water_min_coordinates = water_min_layer.pixel_to_world_coordinates(filter_no_data_value=False)
water_min_coordinates
mintemp_dataframe = pd.DataFrame([water_min_coordinates[0], water_min_coordinates[1]]).T
mintemp_dataframe.columns=['decimallatitude', 'decimallongitude']
water_min_matrix = water_min_reader.read(1)
mintemp_dataframe['MinT'] = water_min_matrix.reshape(np.product(water_min_matrix.shape))
mintemp_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
mintemp_dataframe.head()
mintemp_dataframe.tail()
Explanation: 5.4 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a minimum temperature pixel map
End of explanation
water_max_layer =  ClimateLayer(file_path="../data/watertemp/max_wt_2000.tif") 
water_max_reader = water_max_layer.load_data()
# HERE: should we ignore cells with no-data values for temperature? They are set to a really big negative number
# for now we keep them, otherwise could be NaN
water_max_coordinates = water_max_layer.pixel_to_world_coordinates(filter_no_data_value=False)
maxtemp_dataframe = pd.DataFrame([water_max_coordinates[0], water_max_coordinates[1]]).T
maxtemp_dataframe.columns=['decimallatitude', 'decimallongitude']
water_max_matrix = water_max_reader.read(1)
maxtemp_dataframe['MaxT'] = water_max_matrix.reshape(np.product(water_max_matrix.shape))
maxtemp_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
maxtemp_dataframe.head()
maxtemp_dataframe.tail()
Explanation: 5.5 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a maximum temperature pixel map
End of explanation
water_mean_layer =  ClimateLayer(file_path="../data/watertemp/mean_wt_2000.tif") 
water_mean_reader = water_mean_layer.load_data()
# HERE: should we ignore cells with no-data values for temperature? They are set to a really big negative number
# for now we keep them, otherwise could be NaN
water_mean_coordinates = water_mean_layer.pixel_to_world_coordinates(filter_no_data_value=False)
meantemp_dataframe = pd.DataFrame([water_mean_coordinates[0], water_mean_coordinates[1]]).T
meantemp_dataframe.columns=['decimallatitude', 'decimallongitude']
water_mean_matrix = water_mean_reader.read(1)
meantemp_dataframe['MeanT'] = water_mean_matrix.reshape(np.product(water_mean_matrix.shape))
meantemp_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
meantemp_dataframe.head()
meantemp_dataframe.tail()
# merge base with presences
merged = base_dataframe.combine_first(presences_dataframe)
merged.head()
merged.tail()
# merge based+presences with pseudo-absences
# merged2 = pd.merge(merged1, pseudo_absences_dataframe, on=["decimallatitude", "decimallongitude", salmo_trutta.name_species], how="outer")
merged = merged.combine_first(pseudo_absences_dataframe)
Explanation: 5.6 Get arrays of coordinates (latitude/longitude) for each cell (middle point) in a mean temperature pixel map
End of explanation
merged.head()
merged.tail()
# merge base+presences+pseudo-absences with min temperature
#merged3 = pd.merge(merged2, mintemp_dataframe, on=["decimallatitude", "decimallongitude"], how="outer")
merged = merged.combine_first(mintemp_dataframe)
merged.head()
merged.tail()
# merged4 = pd.merge(merged3, maxtemp_dataframe, on=["decimallatitude", "decimallongitude"], how="outer")
merged = merged.combine_first(maxtemp_dataframe)
merged.head()
merged.tail()
# merged5 = pd.merge(merged4, meantemp_dataframe, on=["decimallatitude", "decimallongitude"], how="outer")
merged = merged.combine_first(meantemp_dataframe)
merged.tail()
merged.to_csv("../data/fish/selection/salmo_trutta.csv")
merged[merged['Salmo trutta']==0].shape[0] # should be equal to number of pseudo absences below
pseudo_absence_coordinates[0].shape[0]
merged[merged['Salmo trutta']==1].shape[0]  # should be equal to number of presences below
presence_coordinates[0].shape[0]
merged[merged['Salmo trutta'].isnull()].shape[0] # all that's left
360 * 720 == merged[merged['Salmo trutta']==0].shape[0] + merged[merged['Salmo trutta']==1].shape[0] + merged[merged['Salmo trutta'].isnull()].shape[0]
# == all pixels in 360 x 720 matrix
Explanation: http://pandas.pydata.org/pandas-docs/stable/merging.html
For this, use the combine_first method.
Note that this method only takes values from the right DataFrame if they are missing in the left DataFrame. A related method, update, alters non-NA values inplace
End of explanation
# Download from Google Drive: https://drive.google.com/open?id=0B9cazFzBtPuCaW0wRkk2N0g5d1k
lepidomeda_mollispinis = IUCNSpecies(name_species='Lepidomeda mollispinis')
lepidomeda_mollispinis.load_shapefile("../data/fish/selection/lepidomeda_mollispinis")
rasterized_lm = lepidomeda_mollispinis.rasterize(raster_file="./lepidomeda_mollispinis_full.tif", pixel_size=0.5)
plt.figure(figsize=(25,20))
plt.imshow(rasterized_lm, cmap="hot", interpolation="none")
selected_layers_lm, pseudo_absences_lm = biomes_adf.sample_pseudo_absences(species_raster_data=rasterized_lm, continents_raster_data=continents_rasters, number_of_pseudopoints=1000)
plt.figure(figsize=(25,20))
plt.imshow(selected_layers_lm, cmap="hot", interpolation="none")
plt.figure(figsize=(25,20))
plt.imshow(pseudo_absences_lm, cmap="hot", interpolation="none")
presence_coordinates_lm = lepidomeda_mollispinis.pixel_to_world_coordinates()
presences_dataframe = pd.DataFrame([presence_coordinates_lm[0], presence_coordinates_lm[1]]).T
presences_dataframe.columns=['decimallatitude', 'decimallongitude']
presences_dataframe[lepidomeda_mollispinis.name_species] = 1 # fill presences with 1's
presences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
presences_dataframe.head()
presences_dataframe.tail()
pseudo_absence_coordinates_lm = biomes_adf.pixel_to_world_coordinates(raster_data=pseudo_absences_lm)
pseudo_absences_dataframe = pd.DataFrame([pseudo_absence_coordinates_lm[0], pseudo_absence_coordinates_lm[1]]).T
pseudo_absences_dataframe.columns=['decimallatitude', 'decimallongitude']
pseudo_absences_dataframe[lepidomeda_mollispinis.name_species] = 0
pseudo_absences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
pseudo_absences_dataframe.head()
pseudo_absences_dataframe.tail()
merged1 = merged.combine_first(presences_dataframe)
merged1.tail()
merged1 = merged1.combine_first(pseudo_absences_dataframe)
merged1.tail()
merged1['Lepidomeda mollispinis'].unique()
merged1[merged1['Lepidomeda mollispinis']==0].shape # pseudo-absences
merged1[merged1['Lepidomeda mollispinis']==1].shape # presences
merged1[merged1['Lepidomeda mollispinis'].isnull()].shape
merged1[merged1['Lepidomeda mollispinis'].isnull()].shape[0] + merged1[merged1['Lepidomeda mollispinis']==1].shape[0] + merged1[merged1['Lepidomeda mollispinis']==0].shape[0]
salmo_trutta.get_data().shape_area.sum()
lepidomeda_mollispinis.get_data().shape_area.sum()
Explanation: 6. Repeat with other species
End of explanation
# Download from Google drive: https://drive.google.com/open?id=0B9cazFzBtPuCamEwWlZxV3lBZmc
esox_lucius = IUCNSpecies(name_species='Esox lucius')
esox_lucius.load_shapefile("../data/fish/selection/esox_lucius/")
rasterized_el = esox_lucius.rasterize(raster_file="./esox_lucius_full.tif", pixel_size=0.5)
plt.figure(figsize=(25,20))
plt.imshow(rasterized_el, cmap="hot", interpolation="none")
selected_layers_el, pseudo_absences_el = biomes_adf.sample_pseudo_absences(species_raster_data=rasterized_el, continents_raster_data=continents_rasters, number_of_pseudopoints=1000)
plt.figure(figsize=(25,20))
plt.imshow(selected_layers_el, cmap="hot", interpolation="none")
Explanation: Third species.... (largest IUCN area, also plenty of occurrences in GBIF)
End of explanation
np.where((continents_rasters[2]+rasterized_el)>1)
Explanation: Hmm, why does it take South Afrika into account?
End of explanation
(continents_rasters[2]+rasterized_el)[108,348]
plt.figure(figsize=(25,20))
plt.imshow(pseudo_absences_el, cmap="hot", interpolation="none")
presence_coordinates_el = esox_lucius.pixel_to_world_coordinates()
presences_dataframe = pd.DataFrame([presence_coordinates_el[0], presence_coordinates_el[1]]).T
presences_dataframe.columns=['decimallatitude', 'decimallongitude']
presences_dataframe[esox_lucius.name_species] = 1 # fill presences with 1's
presences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
presences_dataframe.head()
pseudo_absence_coordinates_el = biomes_adf.pixel_to_world_coordinates(raster_data=pseudo_absences_el)
pseudo_absences_dataframe = pd.DataFrame([pseudo_absence_coordinates_el[0], pseudo_absence_coordinates_el[1]]).T
pseudo_absences_dataframe.columns=['decimallatitude', 'decimallongitude']
pseudo_absences_dataframe[esox_lucius.name_species] = 0 # fill pseudo-absences with 0
pseudo_absences_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
pseudo_absences_dataframe.head()
merged2 = merged1.combine_first(presences_dataframe)
merged2.head()
merged2 = merged2.combine_first(pseudo_absences_dataframe)
merged2.head()
merged2['Esox lucius'].unique()
merged2[merged2['Esox lucius']==0].shape # pseudo-absences
merged2[merged2['Esox lucius']==1].shape # presences
merged2[merged2['Esox lucius'].isnull()].shape
merged2[merged2['Esox lucius'].isnull()].shape[0] + merged2[merged2['Esox lucius']==1].shape[0] + merged2[merged2['Esox lucius']==0].shape[0]
merged2.tail()
# rearange columns (nothing critical)
cols = merged2.columns.values
cols1 = [cols[4], cols[2], cols[3], cols[0], cols[1], cols[5]]
merged2 = merged2[cols1]
merged2.to_csv("../data/fish/selection/dataframe_merged_all_touching_false.csv")
merged2.columns.values
Explanation: This above adds two bands with ones and zeroes: The raster for the species, and the raster for the South-Afrika continent (id=2). Then it finds pixel positions where the "summed band" has value >1. (np.where(...))
Result indicates the pixel x/y positions. Basically they overlap at one pixel position. There, the value is 2.
End of explanation
for idx, band in enumerate(continents_rasters):
    continents_coordinates = biomes_adf.pixel_to_world_coordinates(raster_data=band)
    continent_dataframe = pd.DataFrame([continents_coordinates[0], continents_coordinates[1]]).T
    continent_dataframe.columns=['decimallatitude', 'decimallongitude']
    continent_dataframe['Continent'] = idx
    continent_dataframe.set_index(['decimallatitude', 'decimallongitude'], inplace=True, drop=True)
    continent_dataframe.head()
    merged2 = merged2.combine_first(continent_dataframe)
merged2[merged2.Continent==0].shape
np.count_nonzero(continents_rasters[0]) # good!
merged2[merged2.Continent==2].shape
np.count_nonzero(continents_rasters[2]) # good!
merged2.columns.values
# rearange columns again
cols = merged2.columns.values
cols1 = [cols[5], cols[3], cols[4], cols[1], cols[2], cols[6], cols[0]]
merged2 = merged2[cols1]
merged2.head()
merged2.to_csv("../data/fish/selection/dataframe_merged_with_continents_rasterized_2.csv")
np.where((continents_rasters[2]+rasterized_el)>1)
Explanation: Add continents column
End of explanation
<END_TASK>