In [1]:
from iSDM.species import GBIFSpecies
In [2]:
my_species = GBIFSpecies(name_species="Etheostoma_blennioides")
In [3]:
my_data = my_species.load_data("Etheostoma_blennioides2382397.pkl") # this is a reduced dataset
In [4]:
%matplotlib inline
my_species.polygonize().plot() # crashes if NaN coordinates. osgeo's fault; fixed to ignore them
Out[4]:
In [5]:
my_species.data_full.geometry.plot()
Out[5]:
In [42]:
from shapely.geometry import Point, Polygon
# cut the outliers, i.e., overlay with a polygon
my_species.data_full = my_species.data_full[my_species.data_full.geometry.within(Polygon(((-100,30), (-100, 50), (-70, 50),(-70, 30))))]
In [43]:
my_species.polygonize().plot()
Out[43]:
In [44]:
my_species.polygonize(buffer=0.5, simplify_tolerance=0.05).plot()
Out[44]:
In [45]:
my_species.polygonize(buffer=0.3, simplify_tolerance=0.03).plot() # etc
Out[45]:
In [41]:
filtered_data = my_species.data_full.geometry[my_species.data_full.geometry.within(Polygon(((-100,30), (-100, 50), (-70, 50),(-70, 30))))]
filtered_data.plot()
Out[41]:
In [11]:
my_species.polygonize(buffer=0.3, simplify_tolerance=0.03, with_envelope=True).plot() # with_envelope means pixelized
Out[11]:
In [63]:
type(filtered_data)
Out[63]:
In [66]:
my_species.data_full.shape
Out[66]:
In [67]:
filtered_data.shape
Out[67]:
In [56]:
# cut even more
my_species.polygonize(buffer=0.3, simplify_tolerance=0.03, with_envelope=True).geometry.intersection(Polygon(((-90,40), (-90, 50), (-70, 50),(-70, 40)))).plot()
Out[56]:
In [48]:
Polygon(((-100,30), (-100, 50), (-70, 50),(-70, 30)))
Out[48]:
In [57]:
my_species.polygonize(buffer=0.3, simplify_tolerance=0.03, with_envelope=True).geometry.intersection(Polygon(((-90,40), (-90, 50), (-70, 50),(-70, 40))))
Out[57]:
In [108]:
from geopandas import GeoDataFrame
In [109]:
from shapely.geometry import Point
In [121]:
crs = None # converting GBIF pointrecords to geopandas; maybe a good idea to have it unified
geometry = [Point(xy) for xy in zip(my_data['decimallongitude'], my_data['decimallatitude'])]
In [122]:
geo_df = GeoDataFrame(my_data, crs=crs, geometry=geometry)
In [112]:
%matplotlib inline
geo_df.plot(figsize=(20, 25))
Out[112]:
In [113]:
geo_df.centroid.plot(figsize=(20,25)) # nope, still points
Out[113]:
In [123]:
my_data
Out[123]:
In [126]:
geo_df
Out[126]:
In [127]:
type(my_data)
Out[127]:
In [128]:
my_data['decimallatitude'].head()
Out[128]:
In [129]:
my_data.head() # notice now the geometry column is there, automatically appended
Out[129]:
In [79]:
type(my_data)
Out[79]:
In [130]:
geo_df.shape
Out[130]:
In [131]:
my_data.shape
Out[131]:
In [132]:
my_data.shape == geo_df.shape # hmmm
Out[132]:
In [133]:
geo_df.geometry # geometry in GBIF data is POINT of course, geometry column in IUCN is polygon.
# but can't we make a polygon out of points? Like the outer boundary or something?
Out[133]:
In [134]:
geo_df.geometry.area # right, this is why it won't work, its a series of points, not polygons.
Out[134]:
In [135]:
geo_df.geometry[1]
Out[135]:
In [136]:
my_data.geometry[1].buffer(2) #yay
Out[136]:
In [137]:
geo_df.geometry.buffer(2).plot()
Out[137]:
In [138]:
geo_df.buffer(0.5).plot(figsize=(20,25)) # a buffer around each point record; we could draw pseudo-absences from there
Out[138]:
In [139]:
type(geo_df.geometry[1])
Out[139]:
In [140]:
from iSDM.species import IUCNSpecies
turtles = IUCNSpecies(name_species="Acanthochelys_pallidipectoris")
In [141]:
turtles_df = turtles.load_data("/home/daniela/git/iSDM/notebooks/Acanthochelys pallidipectoris75.pkl")
In [142]:
turtles_df
Out[142]:
In [143]:
turtles_df.geometry
Out[143]:
In [144]:
type(turtles_df.geometry)
Out[144]:
In [145]:
type(geo_df.geometry)
Out[145]:
In [146]:
turtles_df.geometry.bounds
Out[146]:
In [147]:
turtles_df.geometry.length
Out[147]:
In [148]:
geo_df.geometry.bounds.head() # bounds == points, nothing changes. but at least useful to get individual x y coords
Out[148]:
In [150]:
turtles_df.representative_point() # could be useful for sampling?
Out[150]:
In [151]:
turtles_df.geometry.plot()
Out[151]:
In [152]:
geo_df.geometry.representative_point().head()
Out[152]:
In [153]:
geo_df.geometry.head() # representative point of a point, well, still.
Out[153]:
In [154]:
turtles_df.geometry.is_valid
Out[154]:
In [155]:
turtles_df.contains(Point(-61,-26))
Out[155]:
In [156]:
turtles_df.geometry.boundary.plot()
Out[156]:
In [157]:
turtles_df.difference(Point(-61, -26)).plot(figsize=(24,26))
Out[157]:
In [158]:
turtles_df.iat[0,8]
Out[158]:
In [159]:
turtles_df.geometry.convex_hull.plot()
Out[159]:
In [160]:
turtles_df.geometry.envelope.plot()
Out[160]:
In [161]:
turtles_df.simplify(tolerance=0.3).plot() # nice!
Out[161]:
In [162]:
turtles_df.crs # !!!! <-- get the coordinate reference system.
Out[162]:
In [164]:
shapely_point = geo_df.geometry[1]
shapely_point.xy
Out[164]:
In [165]:
shapely_point.buffer(2).area
Out[165]:
In [166]:
# ONE WAY OF CREATING A POLYGON
shapely_buffer = shapely_point.buffer(2.0)
shapely_buffer.simplify(0.2, preserve_topology=False)
Out[166]:
In [167]:
type(shapely_buffer) # already a polygon
Out[167]:
In [170]:
# merge with another point
shapely_buffer.simplify(0.2, preserve_topology=True).union( geo_df.geometry[5].buffer(2).simplify(0.2, preserve_topology=False))
Out[170]:
In [171]:
# merge with another point
another_polygon = shapely_buffer.simplify(0.2, preserve_topology=False).union( geo_df.geometry[10].buffer(2).simplify(0.2, preserve_topology=False))
another_polygon.union( geo_df.geometry[100].buffer(2).simplify(0.2, preserve_topology=True))
Out[171]:
In [172]:
shapely_buffer.exterior.coords.xy
Out[172]:
In [173]:
import shapely.ops
In [229]:
third_polygon = geo_df.geometry[30].buffer(1, resolution=160).simplify(0.1, preserve_topology=False)
In [230]:
third_polygon
Out[230]:
In [191]:
fourth_polygon = geo_df.geometry[50].buffer(2).simplify(0.2, preserve_topology=False)
In [192]:
fourth_polygon
Out[192]:
In [193]:
shapely.ops.cascaded_union([another_polygon, third_polygon, fourth_polygon])
Out[193]:
In [194]:
type(shapely.ops.cascaded_union([another_polygon, third_polygon]))
Out[194]:
In [195]:
multipolyon = shapely.ops.cascaded_union([another_polygon, third_polygon, fourth_polygon])
multipolyon.convex_hull
Out[195]:
In [196]:
third_polygon.envelope # or we can "pixelize" every point (with a buffer around) like this, and then combine the envelopes
Out[196]:
In [197]:
multipolyon_envelope = shapely.ops.cascaded_union([another_polygon.envelope, third_polygon.envelope, fourth_polygon.envelope])
In [200]:
multipolyon_envelope
Out[200]:
In [201]:
multipolyon_envelope.convex_hull
Out[201]:
In [233]:
geo_df.geometry.buffer(1, resolution=160).simplify(0.1, preserve_topology=False)
Out[233]:
In [262]:
# lets do it for all
multi_polygon = shapely.ops.cascaded_union(geo_df.geometry.buffer(1).simplify(0.1, preserve_topology=False).tolist())
In [263]:
multi_polygon
Out[263]:
In [264]:
multi_polygon.convex_hull
Out[264]:
In [265]:
multi_polygon_with_envelope = shapely.ops.cascaded_union(geo_df.geometry.buffer(1).simplify(0.1, preserve_topology=False).envelope.tolist())
In [266]:
multi_polygon_with_envelope
Out[266]:
In [267]:
multi_polygon_with_envelope.convex_hull
Out[267]:
In [268]:
shapely.ops.cascaded_union(geo_df.geometry.buffer(1).envelope.tolist()) # without simplify
Out[268]:
In [271]:
no_simplify = shapely.ops.cascaded_union(geo_df.geometry.buffer(1).tolist()) # without simplify
In [272]:
type(no_simplify)
Out[272]:
In [281]:
no_simplify.convex_hull.boundary
Out[281]:
In [283]:
turtles_df.geometry.plot(figsize=(20,20))
Out[283]:
In [284]:
type(turtles_df.geometry)
Out[284]:
In [285]:
type(no_simplify.convex_hull)
Out[285]:
In [307]:
geo_df_geometry = geo_df.copy(deep=True)
In [440]:
geo_df_geometry = geo_df_geometry.buffer(1).simplify(0.1, preserve_topology=False)
In [308]:
geo_df_geometry['geometry'] = geo_df_geometry['geometry'].apply(lambda x: x.buffer(1).simplify(0.1, preserve_topology=False))
In [309]:
geo_df_geometry.plot()
Out[309]:
In [85]:
# expanding each sample point into its polygon of influence (buffer)
geo_df.buffer(1).simplify(0.1, preserve_topology=False).plot() # no need for lambda like above
In [325]:
shapely.ops.cascaded_union(geo_df_geometry.geometry.tolist()) # this is a multipolygon
Out[325]:
In [326]:
from geopandas import GeoDataFrame, GeoSeries
In [330]:
new_series = GeoSeries(shapely.ops.cascaded_union(geo_df_geometry.geometry.tolist()))
In [331]:
new_series.plot()
Out[331]:
In [333]:
new_series # not a multipolygon?
Out[333]:
In [337]:
type(shapely.ops.cascaded_union(geo_df_geometry.geometry))
Out[337]:
In [338]:
import pandas as pd
In [342]:
len(shapely.ops.cascaded_union(geo_df_geometry.geometry)) # HMMM, each island is a separate polygon
Out[342]:
In [347]:
shapely.ops.cascaded_union(geo_df_geometry.geometry)[4]
Out[347]:
In [348]:
shapely.ops.cascaded_union(geo_df_geometry.geometry)[2]
Out[348]:
In [380]:
my_multipolygon = shapely.ops.cascaded_union(geo_df_geometry.geometry)
In [381]:
type(my_multipolygon)
Out[381]:
In [382]:
len(my_multipolygon)
Out[382]:
In [392]:
again_series = GeoSeries(my_multipolygon)
In [393]:
again_series.shape # all in one polygon??
Out[393]:
In [409]:
again_df = GeoDataFrame( geometry=[pol for pol in my_multipolygon]) # no .tolist for MultiPolygon unfortunatelly
In [410]:
again_df.plot() # YESS THANK YOU. Different islands different polygons
Out[410]:
In [434]:
again_df.geometry.representative_point().plot()
Out[434]:
In [435]:
again_df.geometry.convex_hull.plot()
Out[435]:
In [437]:
again_series.convex_hull.plot()
Out[437]:
In [425]:
geo_df_geometry.geometry
Out[425]:
In [443]:
geo_df.geometry
Out[443]:
In [444]:
geo_df = geo_df.buffer(1).simplify(0.1, preserve_topology=False)
In [445]:
geo_df.geometry
Out[445]:
In [64]:
from geopandas import GeoSeries, GeoDataFrame
ajde = GeoDataFrame.from_file("/home/daniela/git/iSDM/data/urban_areas/ne_50m_urban_areas.shp")
In [63]:
ajde
Out[63]:
In [67]:
from geopandas import GeoDataFrame, GeoSeries
In [68]:
GeoDataFrame.from_file("../data/example.json")
Out[68]:
In [70]:
from osgeo import gdal
dem_file = gdal.Open("/home/daniela/git/iSDM/data/tmax1/tmax1.bil")
In [73]:
num_bands = dem_file.RasterCount
band = dem_file.GetRasterBand(1)
data = band.ReadAsArray()
In [74]:
data
Out[74]:
In [75]:
band.GetNoDataValue()
Out[75]:
In [76]:
data.shape # degrees * 10, timezones? 8640/360=24
Out[76]:
In [81]:
data[data>1].shape
Out[81]:
In [86]:
from geopandas.geoseries import *
In [87]:
from geopandas import GeoDataFrame, GeoSeries
In [88]:
from shapely.geometry import Point, Polygon
In [121]:
p1 = Point(.5,.5)
p2 = Point(.5,1)
p3 = Point(1,1)
p4 = Point(2.1,2.1)
p5 = Point(2,2.5)
p6 = Point(5,5)
In [122]:
points = GeoSeries([p1,p2,p3,p4,p5,p6])
In [127]:
poly = GeoSeries([Polygon([(0,0), (0,2), (2,2), (2,0)]),Polygon([(3,3), (3,6), (6,6), (6,3)])])
In [128]:
points.intersects(poly.unary_union)
Out[128]:
In [129]:
points.plot()
Out[129]:
In [130]:
poly.plot()
Out[130]:
In [120]:
type(poly)
Out[120]:
In [12]:
from osgeo import ogr
from osgeo import gdal
geo = gdal.Open("/home/daniela/git/iSDM/data/tmax1/tmax12.bil")
In [13]:
drv= geo.GetDriver()
In [14]:
drv.GetMetadataItem('DMD_LONGNAME')
Out[14]:
In [15]:
geo.GetLayerCount()
Out[15]:
In [16]:
band = geo.GetRasterBand(1)
In [17]:
type(band)
Out[17]:
In [18]:
band.ReadAsArray().shape
Out[18]:
In [19]:
drv.LongName
Out[19]:
In [20]:
drv.GetMetadataItem("MinX")
In [21]:
band.GetMetadataItem("MinX")
In [22]:
band.DataType
Out[22]:
In [23]:
geotransform = geo.GetGeoTransform()
In [24]:
geotransform
Out[24]:
In [25]:
import numpy as np
data = band.ReadAsArray(0, 0, geo.RasterXSize, geo.RasterYSize).astype(np.float)
In [26]:
band.GetNoDataValue()
Out[26]:
In [27]:
data[data>-9999]
Out[27]:
In [28]:
geo1 = ogr.Open("/home/daniela/git/iSDM/data/tmax1/tmax12.hdr")
In [29]:
type(geo1)
Out[29]:
In [30]:
data[data>-9999].min()/10
Out[30]:
In [31]:
data[data>-9999].max()/10
Out[31]:
In [32]:
band.GetBlockSize()
Out[32]:
In [33]:
band.ComputeRasterMinMax()
Out[33]:
In [34]:
band.GetMinimum() # WTF?
Out[34]:
In [35]:
band.GetMaximum() # WTF?
Out[35]:
In [36]:
band.GetRasterCategoryNames()
In [37]:
gdal.GetDataTypeName(band.DataType)
Out[37]:
In [38]:
band.GetMinimum()
Out[38]:
In [39]:
drv.GetMetadata_List()
Out[39]:
In [40]:
band.ComputeBandStats()
Out[40]:
In [41]:
band.ComputeStatistics(True)
Out[41]:
In [42]:
band.GetColorTable()
In [43]:
band.GetDataset()
Out[43]:
In [44]:
import rasterio
In [45]:
with rasterio.open("/home/daniela/git/iSDM/data/tmax1/tmax12.bil") as src:
print(src.width, src.height)
print(src.crs)
print(src.affine)
print(src.count)
print(src.indexes)
In [46]:
dataset = rasterio.open("/home/daniela/git/iSDM/data/tmax1/tmax12.bil")
In [47]:
dataset.affine
Out[47]:
In [48]:
dataset.colorinterp(1)
Out[48]:
In [49]:
dataset.meta
Out[49]:
In [50]:
dataset.shape
Out[50]:
In [51]:
type(dataset)
Out[51]:
In [52]:
type(dataset.read())
Out[52]:
In [53]:
type(band)
Out[53]:
In [54]:
type(data)
Out[54]:
In [55]:
data.shape
Out[55]:
In [56]:
dataset.read().shape
Out[56]:
In [57]:
dataset.read_band(1).shape
Out[57]:
In [58]:
dataset.bounds # YAY
Out[58]:
In [59]:
from rasterio.plot import show, show_hist
In [60]:
import rasterio.plot # my version is 0.25, they are up to 0.35. Maybe because I'm using python 3.5. Upgraded; works now
In [61]:
dataset.dtypes
Out[61]:
In [62]:
from rasterio.rio.insp import stats
In [127]:
%matplotlib inline
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(18, 15))
ax.imshow(dataset.read(1), cmap="hot")
Out[127]:
In [66]:
dataset.read(1)
Out[66]:
In [78]:
my_data = dataset.read(1)
In [88]:
my_data[my_data>dataset.get_nodatavals()[0]]
Out[88]:
In [83]:
my_data.shape
Out[83]:
In [87]:
dataset.get_nodatavals()[0]
Out[87]:
In [92]:
my_data[2000,5000]/10
Out[92]:
In [95]:
from rasterio.transform import Affine
def transform_from_corner(ulx, uly, dx, dy):
return Affine.translation(ulx, uly)*Affine.scale(dx, -dy)
print(transform_from_corner(bounds[0], bounds[3], 1.0/3600, 1.0/3600).to_gdal())
In [96]:
from rasterio.features import rasterize
from shapely.geometry import Polygon, mapping
# image transform
bounds = (119.52, -21.6, 120.90, -20.5)
transform = transform_from_corner(bounds[0], bounds[3], 1.0/3600, 1.0/3600)
# Make raster image, burn in vector data which lies completely inside the bounding box
poly = Polygon(((120, -21), (120.5, -21), (120.5, -21.2), (120, -21.2)))
output = rasterize([poly], transform=transform, out_shape=(3961, 4969))
print(output)
In [126]:
plt.imshow(output, cmap="hot")
Out[126]:
In [125]:
import rasterio
from rasterio.features import rasterize
from rasterio.transform import IDENTITY
geometry = {
'type': 'Polygon',
'coordinates': [[(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]]}
geometry = Polygon([(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]).buffer(3)
rows = cols = 100
with rasterio.drivers():
result = rasterize([geometry], out_shape=(rows, cols))
with rasterio.open(
"test.tif",
'w',
driver='GTiff',
width=cols,
height=rows,
count=1,
dtype=np.uint8,
nodata=0,
crs={'init': "EPSG:4326"}) as out:
out.write(result.astype(np.uint8), indexes=1)
out.close()
result = None
out = None
In [121]:
Polygon([(2, 2), (2, 4.25), (4.25, 4.25), (4.25, 2), (2, 2)]).buffer(3)
Out[121]:
In [ ]: