<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 1. Collect and filter all observations
Recrods with latitude/longitude
Step1:  How many unique species have occurrence records with latitude/longitude?
Step2:  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.
Step3:  How many unique species HAVE records with latitude/longitude, AND date of event (at least year)
Step4:  How many unique species with latitude/longitude, AND event date after 1990?
Step5:  I guess we should keep only observations of type 'OBSERVATION', 'MACHINE_OBSERVATION' and 'HUMAN_OBSERVATION'?
Step6:  2. Create a bias grid at a resolution of 5arcmin (for the moment it's just counts of observations per pixel)
5 arcmin = 5/60 pixel = 0.083333333 => 'height'
Step7:  Normalize using log10?
Step8:  Standardize?
Step9:  Min-Max scaling?
Step10:  numpy array memory map
Step11:  Can read any tif file as memory map??
<ASSISTANT_TASK:>
Python Code:
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] # total number, before filtering out observations that match our criteria
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")
import gc
gc.collect()
from geopandas import GeoDataFrame
final_observations = GeoDataFrame.from_file("../data/bias_grid/final_observations/")
final_observations.head()
Explanation: I guess we should keep only observations of type 'OBSERVATION', 'MACHINE_OBSERVATION' and 'HUMAN_OBSERVATION'?
End of explanation
x_min, y_min, x_max, y_max = -180, -90, 180, 90
pixel_size = 0.0083333333 # changed from 0.083333333 to 30arcsec
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
bias_grid=np.zeros(shape=(y_res, x_res)).astype('int32')
def increase_pixel_value(row):
    bias_grid[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_value(row))
bias_grid.max()
bias_grid.std()
bias_grid.sum()
# is the sum of the bias grid equal to the total number of observations?
bias_grid.sum() == final_observations.shape[0] 
import gc
gc.collect()
Explanation: 2. Create a bias grid at a resolution of 5arcmin (for the moment it's just counts of observations per pixel)
5 arcmin = 5/60 pixel = 0.083333333 => 'height': 2160, 'width': 4320 for a global map
(180/0.083333333) ~ 2160           
(360/0.083333333) ~ 4320
Generate 2D array and use it as a basis for bias grid.
End of explanation
# bias_grid_plus_1 = bias_grid + 1
bias_grid_log = np.log10(bias_grid + 1)
# bias_grid_log[np.isneginf(bias_grid_log)] = 0 # is this a good idea, setting to 0? log10(0) = infinity otherwise
bias_grid_log.max()
bias_grid_log.min()
bias_grid_log.std()
Explanation: Normalize using log10?
End of explanation
bias_grid_standardized = (bias_grid - bias_grid.mean()) / bias_grid.std()
bias_grid_standardized.max()
bias_grid_standardized.min()
bias_grid_standardized.std()
Explanation: Standardize?
End of explanation
bias_grid_minmax_scale = (bias_grid - bias_grid.min()) / (bias_grid.max() - bias_grid.min())
bias_grid_minmax_scale.max()
bias_grid_minmax_scale.min()
bias_grid_minmax_scale.std()
import matplotlib.pyplot as plt
plt.figure(figsize=(25,20))
plt.imshow(bias_grid_log, cmap="hot", interpolation="none")
import pickle
pickle.dump(bias_grid, open("../data/bias_grid/bias_grid_30arcsec.pkl", "wb"))
np.count_nonzero(bias_grid)
np.product(bias_grid.shape)
41881/9331200
bias_grid.shape
import rasterio
from rasterio.transform import Affine
x_min, y_min, x_max, y_max = -180, -90, 180, 90
pixel_size = 0.083333333
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
crs = {'init': "EPSG:4326"}
transform = Affine.translation(x_min, y_max) * Affine.scale(pixel_size, -pixel_size)
with rasterio.open("../data/bias_grid/bias_grid_minmax_scale.tif", 'w', driver='GTiff', width=x_res, height=y_res,
                   count=1,
                   dtype=np.uint16,
                   nodata=0,
                   transform=transform,
                   crs=crs) as out:
    out.write(bias_grid_minmax_scale.astype(np.uint16), indexes=1)
    out.close()
pixel_size
bias_grid_minmax_scale.std()
bias_grid_log.shape
import pickle
pickle.dump(bias_grid_log, open("../data/bias_grid/bias_grid_log_30arcsec.pkl", "wb"), protocol=4)
pickle.dump(bias_grid_standardized, open("../data/bias_grid/bias_grid_standardized.pkl", "wb"))
pickle.dump(bias_grid_minmax_scale, open("../data/bias_grid/bias_grid_minmax_scale.pkl", "wb"))
Explanation: Min-Max scaling?
End of explanation
# bias_grid_mm=np.zeros(shape=(y_res, x_res)).astype('int32')
bias_grid_mm = np.memmap("../data/bias_grid/bias_grid_mm.dat", dtype='int32', mode='w+', shape=(y_res,x_res))
def increase_pixel_value(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_value(row))
bias_grid_mm.flush()
bias_grid_mm.max()
bias_grid_mm.std()
del bias_grid_mm
gc.collect()
fpr = np.memmap("../data/bias_grid/bias_grid_mm.dat", dtype='int32', mode='r', shape=(y_res, x_res))
fpr.max()
fpr.flags.writeable
Explanation: numpy array memory map
End of explanation
anything = np.memmap("/home/daniela/git/iSDM/data/GLWD/downscaled/original_corrected.tif", dtype='uint8',  mode='r', shape=(y_res,x_res))
np.unique(anything)
type(anything)
gc.collect()
del anything
some_data = np.memmap("../data/bias_grid/some_data.tif", dtype='float64', mode='w+', shape=(y_res,x_res))
type(some_data)
some_data.flags.writeable
isinstance(some_data, np.ndarray) 
some_data.shape
some_data[:,:] = 1
some_data._mmap
some_data
Explanation: Can read any tif file as memory map??
End of explanation
<END_TASK>