In [1]:
from enum import Enum
class Source(Enum):
WORLDCLIM = 1
GLOBE = 2
In [2]:
print("This will be it %s " % list(Source.__members__))
In [3]:
import rasterio
tmax_january = rasterio.open("/home/daniela/git/iSDM/data/tmax1/tmax1.bil")
In [4]:
tmax_january.bounds # it's part of the metadata (header) MinX MinY MaxX MaxY
Out[4]:
In [5]:
tmax_january.get_nodatavals() # also metadata
Out[5]:
In [6]:
tmax_january.crs
Out[6]:
In [7]:
tmax_january.height # NROWS metadata
Out[7]:
In [8]:
tmax_january.lnglat()
Out[8]:
In [9]:
tmax_january.res # pixel size? XDIM YDIM metadata
Out[9]:
In [10]:
tmax_january.meta # awesome!
Out[10]:
In [11]:
tmax_january.meta['transform']
Out[11]:
In [12]:
tmax_january.sample(0.4)
Out[12]:
In [13]:
npixels = tmax_january.width * tmax_january.height
for i in tmax_january.indexes:
band = tmax_january.read(i)
print(i, band[band!=tmax_january.nodata].min(), band.max(), band.sum()/npixels)
In [14]:
npixels
Out[14]:
In [15]:
for key in tmax_january.meta.keys():
print(key,":",tmax_january.meta[key])
In [16]:
print(tmax_january.read())
In [17]:
print(tmax_january.meta)
In [18]:
import json
print(json.dumps(tmax_january.meta, indent=2))
In [19]:
import pprint
In [20]:
pp = pprint.PrettyPrinter(depth=6)
In [21]:
print(pp.pprint(tmax_january.meta))
In [ ]:
In [1]:
In [2]:
In [3]:
temperature_max_january.load_data()
Out[3]:
In [4]:
pp.pformat(tmax_january.meta)
In [18]:
from rasterio.warp import calculate_default_transform, reproject, RESAMPLING
import numpy as np
destination = np.zeros(tmax_january.read(1).shape, np.uint8)
source=rasterio.band(tmax_january, 1)
band = tmax_january.read()
In [19]:
affine, width, height = calculate_default_transform(src_crs=tmax_january.crs,
dst_crs=tmax_january.crs,
width=tmax_january.width,
height=tmax_january.height,
left=tmax_january.bounds.left,
bottom=tmax_january.bounds.bottom,
right=tmax_january.bounds.right,
top=tmax_january.bounds.top,
resolution=5) # CHNGE RESOLUTION
In [122]:
destination = np.zeros((height,width), np.int16)
In [123]:
reproject(source=source,
destination=destination, src_transform=tmax_january.affine,
src_crs=tmax_january.crs,
dst_transform=affine,
dst_crs=tmax_january.crs,
resampling=RESAMPLING.nearest)
In [25]:
destination.shape
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(destination, cmap="hot", interpolation="none")
In [104]:
plt.imshow(band[0], cmap="hot", interpolation="none")
Out[104]:
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import logging
import numpy as np
root = logging.getLogger()
root.addHandler(logging.StreamHandler())
from iSDM.environment import ClimateLayer
temperature_max_january = ClimateLayer(file_path="/home/daniela/git/iSDM/data/tmax1/tmax1.bil")
temperature_max_january.load_dataset()
Out[1]:
In [2]:
temperature_max_january.reproject(destination_file="./here.tif", resolution=10, dst_crs='EPSG:4326')
In [3]:
temperature_max_january.get_data().read()
Out[3]:
In [ ]:
In [4]:
import rasterio
here = rasterio.open("./here.tif")
In [5]:
here.read().shape
Out[5]:
In [6]:
plt.imshow(here.read()[0], cmap="hot", interpolation="none")
Out[6]:
In [7]:
original_data = temperature_max_january.load_dataset()
In [8]:
original_data.shape
Out[8]:
In [9]:
here.read_masks()
Out[9]:
In [10]:
here.meta
Out[10]:
In [11]:
here.bounds
Out[11]:
In [12]:
original_data.read().max() # good, same as /home/daniela/git/iSDM/data/tmax1/tmax1.hdr
Out[12]:
In [13]:
original_data.read().min() # this is the nodata value, we should ignore it
Out[13]:
In [14]:
temperature_mean_january = ClimateLayer(file_path="/home/daniela/git/iSDM/data/tmean/tmean_1/w001001.adf")
temperature_mean_january.load_dataset()
Out[14]:
In [57]:
temperature_mean_january = ClimateLayer(file_path="/home/daniela/git/iSDM/data/tmean/tmean_1/w001001x.adf") # same?
temperature_mean_january.load_dataset()
Out[57]:
In [60]:
temperature_mean_january.reproject(destination_file="./reprojected.tif", resolution=1, driver='EHdr')
In [17]:
source_again.affine
In [15]:
source_again.bounds
In [14]:
import rasterio
climate_data = rasterio.open("/home/daniela/git/iSDM/data/tmax1/tmax1.bil")
In [17]:
climate_data.read()
Out[17]:
In [26]:
climate_data.close()
In [33]:
type(climate_data)
Out[33]:
In [32]:
if not climate_dataa:
print("NO")
In [37]:
temperature_max_january.metadata
Out[37]:
In [38]:
temperature_max_january.metadata.update({'something_else':'value'})
In [39]:
temperature_max_january.metadata
Out[39]:
In [41]:
from iSDM import environment
In [42]:
temperature_max_january.source=environment.Source.WORLDCLIM
In [44]:
temperature_max_january.source.name
Out[44]:
In [45]:
temperature_max_january.metadata.update({'units':'deg C * 10'})
In [46]:
temperature_max_january.metadata
Out[46]:
In [47]:
type(temperature_max_january)
Out[47]:
In [54]:
band = temperature_max_january.get_data().read_band(1)
In [55]:
band
Out[55]:
In [66]:
mask = temperature_max_january.data_full.read_mask(1)
In [68]:
mask.shape
Out[68]:
In [69]:
plt.imshow(mask, cmap="hot", interpolation="none")
Out[69]:
In [70]:
mask
Out[70]:
In [72]:
plt.imshow(band, cmap="hot", interpolation="none")
Out[72]:
In [85]:
band[band>=0]
Out[85]:
In [86]:
mask[mask>=0]
Out[86]:
In [87]:
temperature_max_january.data_full.nodatavals
Out[87]:
In [50]:
# worldclim data
temperature_mean_january = ClimateLayer(file_path="/home/daniela/git/iSDM/data/tmean/tmean_1/w001001.adf")
temperature_mean_january.load_dataset()
Out[50]:
In [51]:
temperature_mean_january.get_data().read(1).shape
Out[51]:
In [56]:
temperature_mean_january.reproject(destination_file="./reprojected.tif", resolution=1, dst_crs='EPSG:4326', driver="AIG")
In [21]:
src = rasterio.open("/home/daniela/git/iSDM/data/tmean/tmean_1/w001001.adf")
In [22]:
src
Out[22]:
In [23]:
from rasterio.warp import calculate_default_transform, RESAMPLING
In [24]:
affine, width, height = calculate_default_transform(src_crs=src.crs,dst_crs=src.crs, width=src.width, height=src.height,left=src.bounds.left,bottom=src.bounds.bottom,right=src.bounds.right,top=src.bounds.top)
In [25]:
affine
Out[25]:
In [26]:
width
Out[26]:
In [27]:
height
Out[27]:
In [49]:
destination_file="./here.tif"
dst = rasterio.open(destination_file, 'w',driver='GTiff', width=width, height=height, count=src.count, dtype=np.int16)
In [29]:
dst
Out[29]:
In [48]:
rasterio.warp.reproject(source=rasterio.band(src, 1),
destination=rasterio.band(dst, 1),
src_transform=src.affine,
src_crs=src.crs,
dst_transform=affine,
dst_crs=src.crs
)
In [46]:
print(affine)
In [126]:
height
Out[126]:
In [47]:
src.indexes
Out[47]:
In [31]:
dst.write(destination_file)
In [32]:
type(destination_file)
Out[32]:
In [57]:
In [66]:
# GLOBE DEM DATA
dem_a = ClimateLayer(file_path="/home/daniela/git/iSDM/data/all10/g10g")
In [67]:
dem_a.load_dataset()
Out[67]:
In [68]:
the_data = dem_a.get_data().read(1)
In [69]:
the_data.shape
Out[69]:
In [70]:
plt.imshow(the_data, cmap="gray", interpolation="none")
Out[70]:
In [71]:
the_data[the_data!=-500].min()
Out[71]:
In [72]:
dem_a.reproject(destination_file="./here.tif", resolution=1, dst_crs='EPSG:4326')
In [76]:
dem_a.get_data().shape
Out[76]:
In [77]:
dem_a.load_dataset(file_path="./here.tif")
Out[77]:
In [78]:
dem_a.get_data().shape # right, updated. Shall I automatically load the new dataset or?
Out[78]:
In [15]:
import rasterio
here = rasterio.open("./here.tif")
In [16]:
here.read().shape
Out[16]:
In [18]:
plt.imshow(here.read()[0], cmap="gray", interpolation="none") # 1 pixel here = 1/0.00833km = 120km
Out[18]:
In [22]:
dem_a.resolution
Out[22]:
In [25]:
# same resolution as GLOBE DEM data, they match in QGIS. But driver is different so reprojection fails
#same as with tmean
worldclim_alt = ClimateLayer(file_path="/home/daniela/git/iSDM/data/alt/alt/w001001.adf")
In [26]:
worldclim_alt.load_dataset()
Out[26]:
In [27]:
worldclim_alt.get_data().read().shape
Out[27]:
In [28]:
worldclim_alt.resolution
Out[28]:
In [39]:
dem_a.reproject(destination_file="./worldclim2.tif", resolution=1, dst_crs='EPSG:4326', driver="AIG")
In [40]:
import rasterio
worldclim2 = rasterio.open("./here.tif")
In [46]:
worldclim2.read(1).shape
Out[46]:
In [42]:
plt.imshow(worldclim2.read()[0], cmap="gray", interpolation="none") # 1 pixel here = 1/0.00833km = 120km
Out[42]:
In [65]:
worldclim_alt.get_data().read().max()
Out[65]:
In [ ]: