CO2 Sequestration

Functions

SetDirectory[NotebookDirectory[]];

getMap[data_]:= Module[{y,x,lat,lon,lonTicks,latTicks,CO2},		
		{y,x} = Dimensions[Reverse@data[[3]]];
		lat[q_] = Ceiling[Fit[Thread[{data[[1]], Range[y]}],{1,q},q]];
		lon[q_] = Ceiling[Fit[Thread[{data[[2]], Range[x]}],{1,q},q]];
		lonTicks = Table[Thread[{Range[x],Floor[data[[2]]+0.05]}][[i]],{i,Append[Range[1,x,300],x]}];
		latTicks = Table[Thread[{Range[y],Reverse@Floor[data[[1]]]+1}][[i]],{i,Append[Range[1,y,100],y]}];
		CO2[lat0_:-90,latf_:89, lon0_:0,lonf_:359] := Reverse@data[[3,lat[lat0];;lat[latf],lon[lon0];;lon[lonf]]];
		
		Return[{CO2,latTicks,lonTicks}]
		]
		
getAcc[y_]:= Module[{data,n,m,ds,dsAcc,p,coords,em,M},
	data = Data[[y]];
	{n,m} = data[[3]]//Dimensions;
	ds = Table[{i,j}->data[[3]][[i,j]],{i,1,n},{j,1,m}];
	ds = Sort[<|Flatten[ds]|>,Greater];
	dsAcc = Accumulate[Values@ds];
	
	Return[{data,ds,dsAcc}]
	]

getTop[data_,ds_,dsAcc_,ρ_]:= Module[{p,coords,em,M},
	{{p}} = Position[dsAcc/Total[Values@ds],Extract[Nearest[dsAcc/Total[Values@ds],ρ],1]];
	coords = Keys[ds[[1;;p]]];
	em = Values[ds[[1;;p]]];
	M = Reverse@SparseArray[Thread[coords-> em],Dimensions[data[[3]]]];
	
	Return[{coords,em,M}]
]

Import Data

years = Range[1970,2012,1];
Data = Table[Import[FileNameJoin[{"CO_Data",ToString[year]<>".nc"}],{"Datasets",{"lat","lon","emi_co"}}],{year,years}];

Generate Maps

maps = Table[
		data = Data[[year]]; 
		{CO2,latTicks,lonTicks} = getMap[data];
		x = CO2[];
		y = 1/(1+Log[x])/.Indeterminate->0.01;
		MatrixPlot[y,ImageSize->Large,FrameTicks->{latTicks,lonTicks},AspectRatio->1/2,ColorFunction->Hue,PlotLabel->years[[year]]],
	{year,Range[1,Length@years]}];

List Animate & GIF

animation = ListAnimate@maps;
Export["CO2.gif",animation]

Get data by percentage

Select Year

{data,ds,dsAcc} = getAcc[Length@years];

Generate data for several percentages

(*Full data*)
{c1,e1,M1} = getTop[data,ds,dsAcc,1];
M1 = 1/(1+Log[Normal@M1])/.Indeterminate->0.01;

(*Smaller percentages*)
{c2,e2,M2} = getTop[data,ds,dsAcc,0.8];
{c3,e3,M3} = getTop[data,ds,dsAcc,0.1];

c_i  are the coordinates of each point
e_i are the emissions of each point
M_i are the matrices of those emissions in their coodinates

Construct maps as graphics

getGraphics[M_]:=Module[{m1,b1,m2,b2,m3,b3},
	m1 = MatrixPlot[M[[1]],ColorFunction->Hue,Frame->False,ImagePadding->0];
	m2 = MatrixPlot[M[[2]],ColorFunction->(Blend[{Transparent,Transparent,Red},#]&),Frame->False,ImagePadding->0];
	m3 = MatrixPlot[M[[3]],ColorFunction->(Blend[{Transparent,Transparent,Black},#]&),Frame->False,ImagePadding->0];
	
	b1 = Graphics@{Texture[m1],Polygon[{ImageScaled/@{{0,0},{1,0},{1,1},{0,1}}},VertexTextureCoordinates->{{0,0},{1,0},{1,1},{0,1}}]};
	b2 = Graphics@{Texture[m2],Polygon[{ImageScaled/@{{0,0},{1,0},{1,1},{0,1}}},VertexTextureCoordinates->{{0,0},{1,0},{1,1},{0,1}}]};
	b3 = Graphics@{Texture[m3],Polygon[{ImageScaled/@{{0,0},{1,0},{1,1},{0,1}}},VertexTextureCoordinates->{{0,0},{1,0},{1,1},{0,1}}]};
	
	Return[{b1,b2,b3}]
	]

{b1,b2,b3} = getGraphics[{M1,M2,M3}];

Show[{b1,b3},AspectRatio->1/2,ImageSize->Large]

Or... with magic

Extract coordinates of the top 10% emissions (variable c3)  from latitude and longitude lists (data[[1]] and data[[2]]) .

coords = Table[{data[[1]][[c3[[i,1]]]],data[[2]][[c3[[i,2]]]]},{i,1,Length@e3}];

Identify the greatest place of CO2 emissions

GeoIdentify["AdministrativeDivision",GeoPosition[coords[[1]]]]

Plot using coordinates

Plot coordinates in world

GeoListPlot[GeoPosition[coords], GeoRange->"World",GeoBackground->"ReliefMap",GeoProjection-> "Equirectangular"]

The first 10 places

GeoListPlot[GeoMarker/@GeoPosition/@coords[[;;10]],GeoRange->"World", GeoBackground->"ContourMap"]

Plot by countries

GeoHistogram[GeoPosition[coords],"Country",GeoRange->"World",PlotLegends->Automatic]

GeoHistogram[GeoPosition[coords],GeoRange->"World",GeoBackground->"CountryBorders",GeoProjection->"WinkelTripel"]

GeoHistogram[GeoPosition[coords],"Rectangle",GeoRange->"World",GeoBackground->"Satellite",GeoProjection->"Albers"]

GeoSmoothHistogram[GeoPosition[coords],GeoRange->"World",GeoBackground->"Coastlines",GeoProjection->"LambertAzimuthal"]

GeoBubbleChart[GeoPosition[coords]->e3,GeoRange->"World",GeoBackground->"Coastlines",GeoProjection->"Mercator"]