In [1]:
import logging
import numpy as np
import pandas as pd
root = logging.getLogger()
root.addHandler(logging.StreamHandler())
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
from iSDM.environment import ContinentsLayer
from iSDM.environment import Source
In [3]:
continents = ContinentsLayer(file_path="../data/continents/continent.shp", source=Source.ARCGIS)
In [4]:
continents.load_data()
In [5]:
continents.get_data()
Out[5]:
In [6]:
fig, ax = plt.subplots(1,1, figsize=(30,20))
continents.data_full.plot(column="continent", colormap="hsv")
Out[6]:
In [7]:
continents_rasters = continents.rasterize(raster_file="../data/continents/continents_raster.tif", pixel_size=0.5, all_touched=True)
In [8]:
continents_rasters.shape # stacked raster with 8 bands, one for each continent.
Out[8]:
In [24]:
plt.figure(figsize=(25,20))
plt.imshow(continents_rasters[4], cmap="hot", interpolation="none")
Out[24]:
In [10]:
plt.figure(figsize=(25,20))
plt.imshow(continents_rasters[4], cmap="hot", interpolation="none")
Out[10]:
In [11]:
raster_reader = continents.load_raster_data()
In [12]:
raster_reader.read(5) # read a particular band
Out[12]:
In [18]:
set(np.unique(continents_rasters))
Out[18]:
In [17]:
continents_rasters.shape[1:]
Out[17]:
In [19]:
isinstance(continents_rasters, np.ndarray)
Out[19]:
In [20]:
set(np.unique(continents_rasters)) == set({0, 1})
Out[20]:
In [21]:
continents_rasters[1]
Out[21]:
In [25]:
# 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")
In [26]:
rasterized = salmo_trutta.rasterize(raster_file="./salmo_trutta_full.tif", pixel_size=0.5, all_touched=True)
In [45]:
plt.figure(figsize=(25,20))
plt.imshow(rasterized, cmap="hot", interpolation="none")
Out[45]:
In [27]:
rasterized.shape
Out[27]:
In [54]:
(continents_rasters * rasterized)[2].max()
Out[54]:
In [60]:
plt.figure(figsize=(25,20))
plt.imshow(continents_rasters[2] + continents_rasters[0], cmap="hot", interpolation="none")
Out[60]:
In [58]:
Out[58]:
In [129]:
overlayed_continents_species = continents_rasters * rasterized
In [130]:
selected_continents = np.zeros_like(rasterized)
In [131]:
for idx, band in enumerate(overlayed_continents_species):
if band.max()==1:
selected_continents +=continents_rasters[idx]
In [121]:
plt.figure(figsize=(25,20))
plt.imshow(selected_continents, cmap="hot", interpolation="none")
Out[121]:
In [118]:
selected_continents.max() # where the continents touch, we have overlap!
Out[118]:
In [122]:
plt.figure(figsize=(25,20))
plt.imshow(selected_continents==2, cmap="hot", interpolation="none")
Out[122]:
In [127]:
selected_continents[selected_continents > 1]=1
In [128]:
plt.figure(figsize=(25,20))
plt.imshow(selected_continents, cmap="hot", interpolation="none")
Out[128]:
In [135]:
if continents_rasters is not None:
print("YES")
In [176]:
continents_rasters = continents_rasters_copy.copy()
In [139]:
continents_rasters_copy = continents_rasters.copy()
In [177]:
continents_rasters[0] = continents_rasters[0] + continents_rasters[2]
In [178]:
continents_rasters = np.delete(continents_rasters, 2, 0) # delete layer 2 (Europe, previously merged with layer 0==Asia)
In [179]:
continents_rasters.shape
Out[179]:
In [182]:
continents_rasters[0].max() # where the continents touch, we have overlap! that's why max is not 1, but 2.
Out[182]:
In [191]:
sum_continents = continents_rasters[0]+continents_rasters[1]+continents_rasters[2]+continents_rasters[3] + continents_rasters[4] + continents_rasters[5]+continents_rasters[6]
In [195]:
plt.figure(figsize=(25,20))
plt.imshow(sum_continents, cmap="hot", interpolation="none")
Out[195]:
In [184]:
continents.get_data().continent
Out[184]:
In [187]:
np.sum(continents_rasters)
Out[187]:
In [ ]: