Rasters¶
Data Exploration¶
Data Source - Forests
Data Source - Surface_Temperature
The data is quite old and shows forests from 2003, which is obviously quite a long time ago, however the methods used here will be applicable to other datasets.
- The data is stored in a tif file and the resolution is 30 arcseconds (~1 square km).
- In each gridpoint of the tif file there is a value between 0–100 corresponding to the percentage tree cover within that 1 square km area.
- Water bodies are given a value of 254 and no data is given 255.
In [1]:
import rasterio
import numpy as np
forests_file = rasterio.open('../resources/gm_ve_v1.tif')
forests = forests_file.read()
print(np.amin(forests))
print(np.amax(forests))
print(len(np.unique(forests)))
0 254 102
In [2]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(facecolor='#FCF6F5FF', figsize=(10,10))
ax.patch.set_facecolor('#FCF6F5FF')
ax.imshow(forests[0],cmap='Greens',interpolation='nearest')
ax.axis('off')
plt.show()
- The obvious issue is that the oceans have been coloured green and there are clearly no forests there.
- The data is scaled between 0–100 for the percentage tree coverage but there is an additional data point for the water bodies — 254.
- When plotting the data matplotlib creates a colourmap with colours scaled linearly between the smallest value (0) and largest value (254) in the dataset.
- The data that we actually care about is mapped to the first 101 colours (~40% of the colourmap), there are 153 colours that are completely unused (101–253) and one colour is used for the oceans (254).
- So ~60% of the colours in the colourmap are not used. This is illustrated visually below.
In [3]:
from matplotlib.colors import ListedColormap, LinearSegmentedColormap, BoundaryNorm
from matplotlib import cm
cmap = cm.get_cmap('Greens', 254)
newcolors = cmap(np.linspace(0, 1, 254))
background_colour = np.array([0.9882352941176471,
0.9647058823529412,
0.9607843137254902,
1.0])
newcolors[:1, :] = background_colour
newcmp_forests = ListedColormap(newcolors)
bounds = np.arange(254)
norm_forests = BoundaryNorm(bounds, newcmp_forests.N)
gradient = np.linspace(0, 1, 254)
gradient = np.vstack((gradient, gradient))
fig, ax = plt.subplots()
ax.imshow(gradient, aspect='auto', cmap=newcmp_forests)
ax.get_yaxis().set_visible(False)
ax.axvline(100, ls="--", c='black', )
plt.text(105, 0.7, "Data Ends", rotation=90)
plt.show()
In [4]:
forests[0][forests[0] == 254] = 0.0
ourcmap = cm.get_cmap('Greens', 101)
newcolors = ourcmap(np.linspace(0, 1, 101))
background_colour = np.array([0.9882352941176471, 0.9647058823529412, 0.9607843137254902, 1.0])
newcolors[:1, :] = background_colour
newcmp_forests = ListedColormap(newcolors)
In [5]:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(facecolor='#FCF6F5FF', figsize=(10,10))
ax.patch.set_facecolor('#FCF6F5FF')
ax.imshow(forests[0],
cmap=newcmp_forests,
interpolation='nearest')
ax.axis('off')
plt.show()
In [6]:
import pyproj
from pyproj import CRS
for x, y in pyproj.pj_list.items():
print(x, y)
adams_hemi Adams Hemisphere in a Square adams_ws1 Adams World in a Square I adams_ws2 Adams World in a Square II aea Albers Equal Area aeqd Azimuthal Equidistant affine Affine transformation airy Airy aitoff Aitoff alsk Modified Stereographic of Alaska apian Apian Globular I august August Epicycloidal axisswap Axis ordering bacon Bacon Globular bertin1953 Bertin 1953 bipc Bipolar conic of western hemisphere boggs Boggs Eumorphic bonne Bonne (Werner lat_1=90) calcofi Cal Coop Ocean Fish Invest Lines/Stations cart Geodetic/cartesian conversions cass Cassini cc Central Cylindrical ccon Central Conic cea Equal Area Cylindrical chamb Chamberlin Trimetric collg Collignon col_urban Colombia Urban comill Compact Miller crast Craster Parabolic (Putnins P4) defmodel Deformation model deformation Kinematic grid shift denoy Denoyer Semi-Elliptical eck1 Eckert I eck2 Eckert II eck3 Eckert III eck4 Eckert IV eck5 Eckert V eck6 Eckert VI eqearth Equal Earth eqc Equidistant Cylindrical (Plate Carree) eqdc Equidistant Conic euler Euler etmerc Extended Transverse Mercator fahey Fahey fouc Foucaut fouc_s Foucaut Sinusoidal gall Gall (Gall Stereographic) geoc Geocentric Latitude geocent Geocentric geogoffset Geographic Offset geos Geostationary Satellite View gins8 Ginsburg VIII (TsNIIGAiK) gn_sinu General Sinusoidal Series gnom Gnomonic goode Goode Homolosine gs48 Modified Stereographic of 48 U.S. gs50 Modified Stereographic of 50 U.S. guyou Guyou hammer Hammer & Eckert-Greifendorff hatano Hatano Asymmetrical Equal Area healpix HEALPix rhealpix rHEALPix helmert 3(6)-, 4(8)- and 7(14)-parameter Helmert shift hgridshift Horizontal grid shift horner Horner polynomial evaluation igh Interrupted Goode Homolosine igh_o Interrupted Goode Homolosine Oceanic View imw_p International Map of the World Polyconic isea Icosahedral Snyder Equal Area kav5 Kavrayskiy V kav7 Kavrayskiy VII krovak Krovak labrd Laborde laea Lambert Azimuthal Equal Area lagrng Lagrange larr Larrivee lask Laskowski lonlat Lat/long (Geodetic) latlon Lat/long (Geodetic alias) latlong Lat/long (Geodetic alias) longlat Lat/long (Geodetic alias) lcc Lambert Conformal Conic lcca Lambert Conformal Conic Alternative leac Lambert Equal Area Conic lee_os Lee Oblated Stereographic loxim Loximuthal lsat Space oblique for LANDSAT mbt_s McBryde-Thomas Flat-Polar Sine (No. 1) mbt_fps McBryde-Thomas Flat-Pole Sine (No. 2) mbtfpp McBride-Thomas Flat-Polar Parabolic mbtfpq McBryde-Thomas Flat-Polar Quartic mbtfps McBryde-Thomas Flat-Polar Sinusoidal merc Mercator mil_os Miller Oblated Stereographic mill Miller Cylindrical misrsom Space oblique for MISR moll Mollweide molobadekas Molodensky-Badekas transformation molodensky Molodensky transform murd1 Murdoch I murd2 Murdoch II murd3 Murdoch III natearth Natural Earth natearth2 Natural Earth 2 nell Nell nell_h Nell-Hammer nicol Nicolosi Globular nsper Near-sided perspective nzmg New Zealand Map Grid noop No operation ob_tran General Oblique Transformation ocea Oblique Cylindrical Equal Area oea Oblated Equal Area omerc Oblique Mercator ortel Ortelius Oval ortho Orthographic pconic Perspective Conic patterson Patterson Cylindrical peirce_q Peirce Quincuncial pipeline Transformation pipeline manager poly Polyconic (American) pop Retrieve coordinate value from pipeline stack push Save coordinate value on pipeline stack putp1 Putnins P1 putp2 Putnins P2 putp3 Putnins P3 putp3p Putnins P3' putp4p Putnins P4' putp5 Putnins P5 putp5p Putnins P5' putp6 Putnins P6 putp6p Putnins P6' qua_aut Quartic Authalic qsc Quadrilateralized Spherical Cube robin Robinson rouss Roussilhe Stereographic rpoly Rectangular Polyconic s2 S2 sch Spherical Cross-track Height set Set coordinate value sinu Sinusoidal (Sanson-Flamsteed) somerc Swiss. Obl. Mercator stere Stereographic sterea Oblique Stereographic Alternative gstmerc Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion) tcc Transverse Central Cylindrical tcea Transverse Cylindrical Equal Area times Times tinshift Triangulation based transformation tissot Tissot tmerc Transverse Mercator tobmerc Tobler-Mercator topocentric Geocentric/Topocentric conversion tpeqd Two Point Equidistant tpers Tilted perspective unitconvert Unit conversion ups Universal Polar Stereographic urm5 Urmaev V urmfps Urmaev Flat-Polar Sinusoidal utm Universal Transverse Mercator (UTM) vandg van der Grinten (I) vandg2 van der Grinten II vandg3 van der Grinten III vandg4 van der Grinten IV vitk1 Vitkovsky I vgridshift Vertical grid shift wag1 Wagner I (Kavrayskiy VI) wag2 Wagner II wag3 Wagner III wag4 Wagner IV wag5 Wagner V wag6 Wagner VI wag7 Wagner VII webmerc Web Mercator / Pseudo Mercator weren Werenskiold I wink1 Winkel I wink2 Winkel II wintri Winkel Tripel xyzgridshift Geocentric grid shift
In [7]:
import rioxarray as rxr
from rasterio.crs import CRS
org_surface_timp_file = rasterio.open('../resources/MOD_LSTD_M_2021-01-01_rgb_3600x1800.TIFF')
org_surface_temp = org_surface_timp_file.read()
surface_temp_file = rxr.open_rasterio('../resources/MOD_LSTD_M_2021-01-01_rgb_3600x1800.TIFF', masked=True).squeeze()
print("Original projection: ", surface_temp_file.rio.crs)
crs_rob = CRS.from_string('+proj=robin')
surface_temp = surface_temp_file.rio.reproject(crs_rob)
print("New projection: ", surface_temp.rio.crs)
surface_temp = surface_temp.to_numpy()
Original projection: EPSG:4326 New projection: PROJCS["unknown",GEOGCS["unknown",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Robinson"],PARAMETER["longitude_of_center",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
Exercise¶
- Plot the new data and the old data to show the difference in projections
- Fix the colourmap issues
- Generate a new projection
In [8]:
fig = plt.figure(facecolor='#FCF6F5FF', figsize=(10,10))
ax1 = plt.subplot(1,2,1)
ax1.patch.set_facecolor('#FCF6F5FF')
ax1.imshow(org_surface_temp[0], cmap='gnuplot', interpolation='nearest')
ax1.axis('off')
ax2 = plt.subplot(1,2,2)
ax2.patch.set_facecolor('#FCF6F5FF')
ax2.imshow(surface_temp, cmap='gnuplot', interpolation='nearest')
ax2.axis('off')
plt.show()
In [9]:
print(np.amin(org_surface_temp), np.amax(org_surface_temp))
print(org_surface_temp[0,0])
org_surface_temp[org_surface_temp == 255] = 0
surface_temp[surface_temp == 255] = 0
cmap = cm.get_cmap('Reds', 255)
newcolors = cmap(np.linspace(0, 1, 255))
bg_color = np.array([0.9882352941176471, 0.9647058823529412, 0.9607843137254902, 1.0])
newcolors[:1, :] = bg_color
newcmp = ListedColormap(newcolors)
bounds = np.arange(256)
norm = BoundaryNorm(bounds, newcmp.N)
0 255 [255 255 255 ... 255 255 255]
In [10]:
fig = plt.figure(facecolor='#FCF6F5FF', figsize=(10,10))
ax1 = plt.subplot(1,2,1)
ax1.patch.set_facecolor('#FCF6F5FF')
ax1.imshow(org_surface_temp[0], cmap=newcmp, norm=norm, interpolation='nearest')
ax1.axis('off')
ax2 = plt.subplot(1,2,2)
ax2.patch.set_facecolor('#FCF6F5FF')
ax2.imshow(surface_temp, cmap=newcmp, norm=norm, interpolation='nearest')
ax2.axis('off')
plt.show()
In [11]:
import rioxarray as rxr
from rasterio.crs import CRS
surface_temp_file = rxr.open_rasterio('../resources/MOD_LSTD_M_2021-01-01_rgb_3600x1800.TIFF', masked=True).squeeze()
print("Original projection: ", surface_temp_file.rio.crs)
crs_rob = CRS.from_string('+proj=merc')
surface_temp = surface_temp_file.rio.reproject(crs_rob)
print("New projection: ", surface_temp.rio.crs)
surface_temp = surface_temp.to_numpy()
Original projection: EPSG:4326 New projection: EPSG:3395
In [12]:
surface_temp[surface_temp == 255] = 0
fig, ax = plt.subplots(facecolor='#FCF6F5FF', figsize=(10,10))
ax.patch.set_facecolor('#FCF6F5FF')
ax.imshow(surface_temp, cmap=newcmp, norm=norm, interpolation='nearest')
ax.axis('off')
ax.set_ylim(2200, 1800)
plt.show()
