ROOT logo

From $ROOTSYS/tutorials/hist/th2polyEurope.C

//This tutorial illustrates how to create an histogram with polygonal
//bins (TH2Poly), fill it and draw it. The initial data are stored
//in TMultiGraphs. They represent the european countries.
//The histogram filling is done according to a Mercator projection,
//therefore the bin contains should be proportional to the real surface
//of the countries.
//
//The initial data have been downloaded from: http://www.maproom.psu.edu/dcw/
//This database was developed in 1991/1992 and national boundaries reflect
//political reality as of that time.
//
//The script is shooting npoints (script argument) randomly over the Europe area.
//The number of points inside the countries should be proportional to the country surface
//The estimated surface is compared to the surfaces taken from wikipedia.
//Author: Olivier Couet

void th2polyEurope(Int_t npoints=500000)
{
   Int_t i,j;
   Double_t lon1 = -25;
   Double_t lon2 =  35;
   Double_t lat1 =  34;
   Double_t lat2 =  72;
   Double_t R = (lat2-lat1)/(lon2-lon1);
   Int_t W    = 800;
   Int_t H    = (Int_t)(R*800);
   gStyle->SetTitleX(0.2);
   gStyle->SetStatX(0.28);
   gStyle->SetStatY(0.45);
   gStyle->SetStatW(0.15);

   // Canvas used to draw TH2Poly (the map)
   TCanvas *ce = new TCanvas("ce", "ce",0,0,W,H);
   ce->ToggleEventStatus();
   ce->SetGridx();
   ce->SetGridy();

   // Real surfaces taken from Wikipedia.
   const Int_t nx = 36;
   // see http://en.wikipedia.org/wiki/Area_and_population_of_European_countries
   char *countries[nx] = { "france",     "spain",  "sweden",  "germany",       "finland",
                           "norway",     "poland", "italy",   "yugoslavia",    "united_kingdom", 
			   "romania",    "belarus","greece",  "czechoslovakia","bulgaria", 
			   "iceland",    "hungary","portugal","austria",       "ireland", 
			   "lithuania",  "latvia", "estonia", "denmark",       "netherlands",
                           "switzerland","moldova","belgium", "albania",       "cyprus", 
			   "luxembourg", "andorra","malta",   "liechtenstein", "san_marino", 
			   "monaco" };
   Float_t surfaces[nx] = { 547030,        505580,   449964,      357021,        338145, 
                            324220,        312685,   301230,      255438,        244820, 
			    237500,        207600,   131940,      127711,        110910, 
			    103000,         93030,    89242,       83870,         70280, 
			     65200,         64589,    45226,       43094,         41526, 
			     41290,         33843,    30528,       28748,          9250, 
			      2586,           468,      316,         160,            61,
                                 2};

   TH1F *h = new TH1F("h","Countries surfaces (in km^{2})",3,0,3);
   for (i=0; i<nx; i++) h->Fill(countries[i], surfaces[i]);
   h->LabelsDeflate();

   TFile::SetCacheFileDir(".");
   TFile *f;
   f = TFile::Open("http://root.cern.ch/files/europe.root","cacheread");

   TH2Poly *p = new TH2Poly(
             "Europe",
             "Europe (bin contents are normalized to the surfaces in km^{2})",
             lon1,lon2,lat1,lat2);
   p->GetXaxis()->SetNdivisions(520);
   p->GetXaxis()->SetTitle("longitude");
   p->GetYaxis()->SetTitle("latitude");

   p->SetContour(100);

   TMultiGraph *mg;
   TKey *key;
   TIter nextkey(gDirectory->GetListOfKeys());
   while (key = (TKey*)nextkey()) {
      obj = key->ReadObj();
      if (obj->InheritsFrom("TMultiGraph")) {
         mg = (TMultiGraph*)obj;
         p->AddBin(mg);
      }
   }

   TRandom r;
   Double_t longitude, latitude;
   Double_t x, y, pi4 = TMath::Pi()/4, alpha = TMath::Pi()/360;
   
   gBenchmark->Start("Partitioning");
   p->ChangePartition(100, 100);
   gBenchmark->Show("Partitioning");

   // Fill TH2Poly according to a Mercator projection.
   gBenchmark->Start("Filling");
   for (i=0; i<npoints; i++) {
      longitude = r.Uniform(lon1,lon2);
      latitude  = r.Uniform(lat1,lat2);
      x         = longitude;
      y         = 38*TMath::Log(TMath::Tan(pi4+alpha*latitude));
      p->Fill(x,y);
   }
   gBenchmark->Show("Filling");

   Int_t nbins = p->GetNumberOfBins();
   Double_t maximum = p->GetMaximum();


   // h2 contains the surfaces computed from TH2Poly.
   TH1F *h2 = h->Clone("h2");
   h2->Reset();
   for (j=0; j<nx; j++) {
      for (i=0; i<nbins; i++) {
         if (strstr(countries[j],p->GetBinName(i+1))) {
            h2->Fill(countries[j],p->GetBinContent(i+1));
            h2->SetBinError(j, p->GetBinError(i+1));
         }
      }
   }

   // Normalize the TH2Poly bin contents to the real surfaces.
   Double_t scale = surfaces[0]/maximum;
   for (i=0; i<nbins; i++) p->SetBinContent(i+1, scale*p->GetBinContent(i+1));

   gStyle->SetOptStat(1111);
   gStyle->SetPalette(1);
   p->Draw("COLZ");

   TCanvas *c1 = new TCanvas("c1", "c1",W+10,0,W-20,H);
   c1->SetRightMargin(0.047);

   Double_t scale = h->GetMaximum()/h2->GetMaximum();

   h->SetStats(0);
   h->SetLineColor(kRed-3);
   h->SetLineWidth(2);
   h->SetMarkerStyle(20);
   h->SetMarkerColor(kBlue);
   h->SetMarkerSize(0.8);
   h->Draw("LP");
   h->GetXaxis()->SetLabelFont(42);
   h->GetXaxis()->SetLabelSize(0.03);
   h->GetYaxis()->SetLabelFont(42);

   h2->Scale(scale);
   Double_t scale2=TMath::Sqrt(scale);
   for (i=0; i<nx; i++) h2->SetBinError(i+1, scale2*h2->GetBinError(i+1));
   h2->Draw("E SAME");
   h2->SetMarkerStyle(20);
   h2->SetMarkerSize(0.8);

   TLegend *leg = new TLegend(0.5,0.67,0.92,0.8,NULL,"NDC");
   leg->SetTextFont(42);
   leg->SetTextSize(0.025);
   leg->AddEntry(h,"Real countries surfaces from Wikipedia (in km^{2})","lp");
   leg->AddEntry(h2,"Countries surfaces from TH2Poly (with errors)","lp");
   leg->Draw();
   leg->Draw();
   
   Double_t wikiSum = h->Integral();
   Double_t polySum = h2->Integral();
   Double_t error = TMath::Abs(wikiSum-polySum)/wikiSum;
   printf("THPoly Europe surface estimation error wrt wikipedia = %f per cent when using %d points\n",100*error,npoints);
}
 th2polyEurope.C:1
 th2polyEurope.C:2
 th2polyEurope.C:3
 th2polyEurope.C:4
 th2polyEurope.C:5
 th2polyEurope.C:6
 th2polyEurope.C:7
 th2polyEurope.C:8
 th2polyEurope.C:9
 th2polyEurope.C:10
 th2polyEurope.C:11
 th2polyEurope.C:12
 th2polyEurope.C:13
 th2polyEurope.C:14
 th2polyEurope.C:15
 th2polyEurope.C:16
 th2polyEurope.C:17
 th2polyEurope.C:18
 th2polyEurope.C:19
 th2polyEurope.C:20
 th2polyEurope.C:21
 th2polyEurope.C:22
 th2polyEurope.C:23
 th2polyEurope.C:24
 th2polyEurope.C:25
 th2polyEurope.C:26
 th2polyEurope.C:27
 th2polyEurope.C:28
 th2polyEurope.C:29
 th2polyEurope.C:30
 th2polyEurope.C:31
 th2polyEurope.C:32
 th2polyEurope.C:33
 th2polyEurope.C:34
 th2polyEurope.C:35
 th2polyEurope.C:36
 th2polyEurope.C:37
 th2polyEurope.C:38
 th2polyEurope.C:39
 th2polyEurope.C:40
 th2polyEurope.C:41
 th2polyEurope.C:42
 th2polyEurope.C:43
 th2polyEurope.C:44
 th2polyEurope.C:45
 th2polyEurope.C:46
 th2polyEurope.C:47
 th2polyEurope.C:48
 th2polyEurope.C:49
 th2polyEurope.C:50
 th2polyEurope.C:51
 th2polyEurope.C:52
 th2polyEurope.C:53
 th2polyEurope.C:54
 th2polyEurope.C:55
 th2polyEurope.C:56
 th2polyEurope.C:57
 th2polyEurope.C:58
 th2polyEurope.C:59
 th2polyEurope.C:60
 th2polyEurope.C:61
 th2polyEurope.C:62
 th2polyEurope.C:63
 th2polyEurope.C:64
 th2polyEurope.C:65
 th2polyEurope.C:66
 th2polyEurope.C:67
 th2polyEurope.C:68
 th2polyEurope.C:69
 th2polyEurope.C:70
 th2polyEurope.C:71
 th2polyEurope.C:72
 th2polyEurope.C:73
 th2polyEurope.C:74
 th2polyEurope.C:75
 th2polyEurope.C:76
 th2polyEurope.C:77
 th2polyEurope.C:78
 th2polyEurope.C:79
 th2polyEurope.C:80
 th2polyEurope.C:81
 th2polyEurope.C:82
 th2polyEurope.C:83
 th2polyEurope.C:84
 th2polyEurope.C:85
 th2polyEurope.C:86
 th2polyEurope.C:87
 th2polyEurope.C:88
 th2polyEurope.C:89
 th2polyEurope.C:90
 th2polyEurope.C:91
 th2polyEurope.C:92
 th2polyEurope.C:93
 th2polyEurope.C:94
 th2polyEurope.C:95
 th2polyEurope.C:96
 th2polyEurope.C:97
 th2polyEurope.C:98
 th2polyEurope.C:99
 th2polyEurope.C:100
 th2polyEurope.C:101
 th2polyEurope.C:102
 th2polyEurope.C:103
 th2polyEurope.C:104
 th2polyEurope.C:105
 th2polyEurope.C:106
 th2polyEurope.C:107
 th2polyEurope.C:108
 th2polyEurope.C:109
 th2polyEurope.C:110
 th2polyEurope.C:111
 th2polyEurope.C:112
 th2polyEurope.C:113
 th2polyEurope.C:114
 th2polyEurope.C:115
 th2polyEurope.C:116
 th2polyEurope.C:117
 th2polyEurope.C:118
 th2polyEurope.C:119
 th2polyEurope.C:120
 th2polyEurope.C:121
 th2polyEurope.C:122
 th2polyEurope.C:123
 th2polyEurope.C:124
 th2polyEurope.C:125
 th2polyEurope.C:126
 th2polyEurope.C:127
 th2polyEurope.C:128
 th2polyEurope.C:129
 th2polyEurope.C:130
 th2polyEurope.C:131
 th2polyEurope.C:132
 th2polyEurope.C:133
 th2polyEurope.C:134
 th2polyEurope.C:135
 th2polyEurope.C:136
 th2polyEurope.C:137
 th2polyEurope.C:138
 th2polyEurope.C:139
 th2polyEurope.C:140
 th2polyEurope.C:141
 th2polyEurope.C:142
 th2polyEurope.C:143
 th2polyEurope.C:144
 th2polyEurope.C:145
 th2polyEurope.C:146
 th2polyEurope.C:147
 th2polyEurope.C:148
 th2polyEurope.C:149
 th2polyEurope.C:150
 th2polyEurope.C:151
 th2polyEurope.C:152
 th2polyEurope.C:153
 th2polyEurope.C:154
 th2polyEurope.C:155
 th2polyEurope.C:156
 th2polyEurope.C:157
 th2polyEurope.C:158
 th2polyEurope.C:159
 th2polyEurope.C:160
 th2polyEurope.C:161
 th2polyEurope.C:162
 th2polyEurope.C:163
 th2polyEurope.C:164
 th2polyEurope.C:165
 th2polyEurope.C:166
thumb
thumb