In [1]:
import logging
root = logging.getLogger()
root.addHandler(logging.StreamHandler())
import matplotlib.pyplot as plt
%matplotlib inline
Loading the data prints also some useful information on the number of raster bands, coordinate reference system, driver, affine transformation, height/width, and bounding box. The availability of this metadata depends on the content of the header file associated with the raster data.
In [2]:
from iSDM.environment import ClimateLayer
worldclim_max_june = ClimateLayer(file_path="../data/tmax1/tmax6.bil")
worldclim_max_june.load_data()
Out[2]:
We can set the source of this data (as additional metadata)
In [3]:
from iSDM.environment import Source
worldclim_max_june.set_source(Source.WORLDCLIM)
This can also be done at the same time while creating the ClimateLayer. Always load the data before reading it
In [4]:
worldclim_max_june = ClimateLayer(file_path="../data/tmax1/tmax6.bil", source=Source.WORLDCLIM)
worldclim_max_june.load_data()
Out[4]:
In [5]:
worldclim_data = worldclim_max_june.read() # read all bands
In [6]:
worldclim_data
Out[6]:
In [7]:
worldclim_data.shape # 3600 rows, 8640 columns, 1 band
Out[7]:
In [8]:
worldclim_data = worldclim_max_june.read(1) # read band 1
In [9]:
worldclim_data.shape # now it's flat matrix datastructure, for that particular band
Out[9]:
In [10]:
worldclim_data
Out[10]:
In [11]:
worldclim_max_june.metadata['nodata'] # this returns the "nodata" value used for cells with no data.
Out[11]:
Close the dataset to release memory, once you are done
In [12]:
worldclim_max_june.close_dataset()
Trying to read a closed dataset ...
In [13]:
worldclim_max_june.read(1) # read band 1
The worldclim_max_june Climate layer instance remembers the file from which it loaded data last time, so you can simply load it again.
In [14]:
worldclim_max_june.load_data()
Out[14]:
... unless of course you decide to change the file from which to load:
In [15]:
worldclim_max_june.load_data("../data/tmax1/tmax1.bil") # let's take january instead of june
Out[15]:
In [16]:
worldclim_max_june.resolution # some more information...
Out[16]:
Python even allows you to freely add properties of this layer object:
In [17]:
worldclim_max_june.mynote = "Random metadata note"
In [18]:
worldclim_max_june.mynote
Out[18]:
In [19]:
worldclim_data = worldclim_max_june.read(1) # read the data again, since we loaded a new dataset (january)
Going back to the band data we loaded, let's find the minimum and maximum values:
In [20]:
worldclim_data.max() # Units = deg C * 10
Out[20]:
In [21]:
worldclim_data.min()
Out[21]:
In [22]:
worldclim_data[worldclim_data!=-9999].min()
Out[22]:
...which is equivalent to:
In [23]:
worldclim_data[worldclim_data!=worldclim_max_june.metadata['nodata']].min()
Out[23]:
Plotting the raster data is just a bit more involved:
In [24]:
fig, ax = plt.subplots(figsize=(15, 15))
value_min = worldclim_data[worldclim_data!=-9999].min()
value_max = worldclim_data.max()
ax.imshow(worldclim_data, cmap="coolwarm", vmin=value_min, vmax=value_max)
Out[24]:
Remember, we only stored the actual pixel values in worldclim_data. It is just a numpy array.
In [25]:
type(worldclim_data)
Out[25]:
In [26]:
worldclim_max_june.reproject(destination_file="../data/tmax1/tmax1_lower_res.bil", resolution=3.0)
The reprojected data is stored in a separate destination_file in order to not overwrite the original data (btw do we want that?). Let's load that data now. Notice the resolution and the height/width
In [27]:
worldclim_max_june.load_data("../data/tmax1/tmax1_lower_res.bil")
Out[27]:
In [28]:
worldclim_data_low_res = worldclim_max_june.read(1)
In [29]:
worldclim_data_low_res.shape # height, width
Out[29]:
Plot the data again. Notice it is "pixelized" i.e., not that smooth because of lower resolution.
In [30]:
fig, ax = plt.subplots(figsize=(15, 15))
value_min = worldclim_data_low_res[worldclim_data_low_res!=-9999].min()
value_max = worldclim_data_low_res.max()
ax.imshow(worldclim_data_low_res, cmap="coolwarm", interpolation="none", vmin=value_min, vmax=value_max)
Out[30]:
In [31]:
from rasterio.warp import RESAMPLING
worldclim_max_june.load_data("../data/tmax1/tmax1.bil")
worldclim_max_june.reproject(destination_file="../data/tmax1/tmax1_lower_res.bil", resolution=3.0, resampling=RESAMPLING.cubic_spline)
Load the data again, and visualize it
In [32]:
worldclim_max_june.load_data("../data/tmax1/tmax1_lower_res.bil")
worldclim_data_low_res = worldclim_max_june.read(1)
In [33]:
fig, ax = plt.subplots(figsize=(15, 15))
value_min = worldclim_data_low_res[worldclim_data_low_res!=-9999].min()
value_max = worldclim_data_low_res.max()
ax.imshow(worldclim_data_low_res, cmap="coolwarm", interpolation="none", vmin=value_min, vmax=value_max)
Out[33]:
The original coordinate reference system is 'epsg:4326'. The target coordinate system may be any of the usual GDAL/OGR forms, complete WKT, PROJ.4, EPSG:n. For more complex examples of reprojection based on new bounds, dimensions, and resolution, using rasterio: https://github.com/mapbox/rasterio/blob/master/docs/reproject.rst and using gdal: http://www.geos.ed.ac.uk/~smudd/TopoTutorials/html/tutorial_raster_conversion.html (rasterio works on top of GDAL)
as well as the command line interface: https://github.com/mapbox/rasterio/blob/master/docs/cli.rst#warp
In [34]:
from iSDM.environment import DEMLayer
dem_layer = DEMLayer(file_path="../data/alt_5m_bil/alt.bil") # altitude data from http://www.worldclim.org/current
dem_layer.load_data()
Out[34]:
In [35]:
from iSDM.environment import Source
dem_layer.set_source(Source.WORLDCLIM)
In [36]:
dem_data = dem_layer.read() # read all bands
dem_data = dem_layer.read(1) # read band 1
In [37]:
dem_data.shape # == height x width
Out[37]:
In [38]:
dem_layer.metadata['nodata']
Out[38]:
In [39]:
dem_data.max() # Units = meters
Out[39]:
In [40]:
dem_data[dem_data!=-9999].min()
Out[40]:
Plot the raster (always calculate the vmin/vmax)
In [41]:
fig, ax = plt.subplots(figsize=(15, 15))
value_min = dem_data[dem_data!=-9999].min()
value_max = dem_data.max()
ax.imshow(dem_data, cmap="seismic", vmin=value_min, vmax=value_max) # or cmap="terrain"?
Out[41]:
Reproject to a lower resolution (the value is the pixel size, in terms of the Units - degrees in this case)
In [42]:
dem_layer.reproject(destination_file="../data/alt_5m_bil/downscaled.bil", resolution=2.0)
Load the projection (shall we make an "overwrite" option to avoid loading every time) to visualize it
In [43]:
dem_layer.load_data("../data/alt_5m_bil/downscaled.bil")
dem_low_res = dem_layer.read(1)
In [44]:
fig, ax = plt.subplots(figsize=(15, 15))
value_min = dem_low_res[dem_low_res!=-9999].min()
value_max = dem_low_res.max()
ax.imshow(dem_low_res, cmap="seismic", interpolation="none", vmin=value_min, vmax=value_max)
Out[44]:
Load back the original data and crop it with a bounding box
In [45]:
dem_layer.load_data("../data/alt_5m_bil/alt.bil")
dem_layer.reproject(destination_file="../data/alt_5m_bil/cropped.bil", left=-50, bottom=42,right=-76, top=54)
Load the cropped data to visualize (notice the BoundingBox values of the loaded data, they match with the values we specified above)
In [46]:
dem_layer.load_data("../data/alt_5m_bil/cropped.bil")
dem_cropped = dem_layer.read(1)
fig, ax = plt.subplots(figsize=(15, 15))
value_min = dem_cropped[dem_cropped!=-9999].min()
value_max = dem_cropped.max()
ax.imshow(dem_cropped, cmap="seismic", interpolation="none", vmin=value_min, vmax=value_max)
Out[46]:
Lower the resolution of the cropped area
In [47]:
dem_layer.reproject(destination_file="../data/alt_5m_bil/cropped_downscaled.bil", resolution=1.0)
In [48]:
dem_layer.load_data("../data/alt_5m_bil/cropped_downscaled.bil")
dem_cropped_downscaled = dem_layer.read(1)
fig, ax = plt.subplots(figsize=(15, 15))
value_min = dem_cropped_downscaled[dem_cropped_downscaled!=-9999].min()
value_max = dem_cropped_downscaled.max()
ax.imshow(dem_cropped_downscaled, cmap="seismic", interpolation="none", vmin=value_min, vmax=value_max)
Out[48]:
In [ ]: