result_with_lat_long = result_with_lat_long[result_with_lat_long.decimallatitude.notnull() & result_with_lat_long.decimallongitude.notnull()]

"""
Explanation: 1. Collect and filter all observations
Recrods with latitude/longitude
End of explanation
"""


result_with_lat_long['species'].unique().size

"""
Explanation: How many unique species have occurrence records with latitude/longitude?
End of explanation
"""


grouped_lat_long_year_or_eventdate = pd.DataFrame()
grouped_lat_long_year_or_eventdate['count'] = result_with_lat_long[result_with_lat_long.eventdate.notnull() | result_with_lat_long.year.notnull()].groupby(['species']).apply(lambda x: x['species'].count())
grouped_lat_long_year_or_eventdate.head(10) # peak at the top 10 only

"""
Explanation: Best to take into account all observations which have either "year" or "eventdate" present. (or both) Let's group them by species name, and count the number of observation records.
End of explanation
"""


result_with_lat_long['species'].unique().size

"""
Explanation: How many unique species HAVE records with latitude/longitude, AND date of event (at least year)
End of explanation
"""


year_or_eventdate_1990 = result_with_lat_long[['species', 'year', 'eventdate', 'basisofrecord', 'decimallatitude', 'decimallongitude']][(result_with_lat_long.year>1990) | (result_with_lat_long.eventdate>"1990")]

grouped_year_or_eventdate_1990 = pd.DataFrame()
grouped_year_or_eventdate_1990['numobservations'] = year_or_eventdate_1990.groupby(['species']).apply(lambda x: x['species'].count())
grouped_year_or_eventdate_1990.shape[0]

year_or_eventdate_1990.basisofrecord.unique()

"""
Explanation: How many unique species with latitude/longitude, AND event date after 1990?
End of explanation
"""


final_selection = year_or_eventdate_1990[(year_or_eventdate_1990.basisofrecord=='OBSERVATION') | (year_or_eventdate_1990.basisofrecord=='HUMAN_OBSERVATION') | (year_or_eventdate_1990.basisofrecord=='MACHINE_OBSERVATION')]

final_selection.species.unique().shape

final_selection

from iSDM.species import GBIFSpecies

all_species = GBIFSpecies(name_species='All')

all_species.set_data(final_selection)

all_species.get_data().species.unique().shape # these many different species

all_species.get_data().shape[0] # 1939675? this many observations satisfying our criteria (after 1990, with the correct observation type)

year_or_eventdate_1990.shape[0] # before filtering out types of observations

all_species.geometrize()

all_species.get_data().species.unique().shape

final_observations = all_species.get_data()[['species', 'year','eventdate', 'basisofrecord','geometry']]

# final_observations.to_file("../data/bias_grid/final_observations", driver="ESRI Shapefile")

"""
Explanation: I guess we should keep only observations of type 'OBSERVATION', 'MACHINE_OBSERVATION' and 'HUMAN_OBSERVATION'?
End of explanation
"""


from geopandas import GeoDataFrame
final_observations = GeoDataFrame.from_file("../data/bias_grid/final_observations/")

final_observations.head()

"""
Explanation: Read all filtered data.
To this point all the operations were to filter out data and finally store the remaining observations in a file.
We can start with the stored file from here.
End of explanation
"""


x_min, y_min, x_max, y_max = -180, -90, 180, 90
pixel_size = 0.083333333 # 0.0083333333
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)

x_res

y_res

"""
Explanation: 2. Create a background grid at a resolution of 30arcsec
30 arcsec = 0.5/60 pixel = 0.0083333333 'height': 21600, 'width': 43200 for a global map
Generate 2D array and use it as a basis for bias grid.
End of explanation
"""


bias_grid_mm = np.memmap("../data/bias_grid/bias_grid_mm_5arcmin.dat", dtype='int32', mode='w+', shape=(y_res,x_res))

def increase_pixel(row):
    bias_grid_mm[np.abs(int((row.y - 90) / pixel_size)),
              np.abs(int((row.x + 180) / pixel_size))]+=1

here = final_observations.geometry.apply(lambda row: increase_pixel(row)) # keeps the memory usage stable

bias_grid_mm.max()

bias_grid_mm.std() # this still eats memory (uses intermediate datastructures, of course)

bias_grid_mm.sum()

# check if the number of all observations is equal to the bias_grid sum of observations
bias_grid_mm.sum() == final_observations.shape[0]

bias_grid_mm.flush()

del bias_grid_mm

# to read it, the following is used:

bias_grid_mm = np.memmap("../data/bias_grid/bias_grid_mm.dat", dtype='int32', mode='r', shape=(y_res,x_res))

bias_grid_mm.sum()

import gc
gc.collect()

"""
Explanation: We will use a memory-mapped biasgrid file to prevent the RAM from exploding. At least while results are computed.
End of explanation
"""


import rasterio
from rasterio.transform import Affine
transform = Affine.translation(x_min, y_max) * Affine.scale(pixel_size, -pixel_size)
crs = {'init': "EPSG:4326"}

transform

with rasterio.open("../data/bias_grid/bias_grid_5arcmin.tif", 'w', driver='GTiff', width=x_res, height=y_res,
                   count=1,
                   dtype=np.uint32,
                   nodata=0,
                   transform=transform,
                   crs=crs) as out:
    out.write(bias_grid_mm.astype(np.uint32), indexes=1)
    out.close()

bias_grid_mm.shape

"""
Explanation: 3. Store as a TIF raster.
End of explanation
"""


