Spmap issue

[Ema Davies]

Hello,

I am working with EU-level data that include different NUTS levels. Is there a way to use spmap in this case? Below is an example of the data (NUTS_RG_03M_2016_3035_LEVL_.shp files are obtained from https://gisco-services.ec.europa.eu/...nuts-2016.html)

* Example generated by -dataex-. For more info, type help dataex
clear
input str5 NUTS_ID float(total x_cc y_cc)
"AT"    2  4631152.5   2728758
"BE1"   1    3924979   3095776
"BE3"   2    3965977   3034637
"BG"    6    5561655 2310162.5
"CY"    0    6426424   1653395
"CZ"    6  4705460.5   2973015
"DE1"   2    4250334   2825138
"DE2"   2  4425542.5   2870462
"DE3"   2  4551963.5 3271179.5
"DE4"   2  4552249.5   3267538
"DE5"   0    4236640   3345796
"DE6"   2  4321787.5   3382180
"DE7"   2  4252034.5   3054615
"DE8"   2    4488061 3408451.5
"DE9"   2    4264924   3295466
"DEA"   0    4151113   3154747
"DEB"   1  4138005.5   2980828
"DEC"   0  4099770.5   2923665
"DED"   2  4555406.5   3109965
"DEE"   2    4437922   3212359
"DEF"   1    4309255   3452599
"DEG"   1    4393193   3088532
"DK"    2  4324183.5   3650497
"EE"    2    5216204   4050577
"EL"    0    5446546   1866405
"EL30"  6    5519687 1747141.6
"EL41"  0    5718016   1894183
"EL42"  0    5781137 1675096.3
"EL43"  1    5677244   1484150
"EL51"  0    5587258 2134529.8
"EL52"  2    5414566   2054253
"EL53"  0    5297153 1992473.8
"EL54"  1    5248981 1900077.8
"EL61"  0    5372253   1911631
"EL62"  0    5232647   1788519
"EL63"  0    5336454   1765899
"EL64"  0    5439385 1824763.5
"EL65"  0    5428266 1675550.5
"ES"    2    3155919   2006099
"ES11"  0  2866097.5 2357738.8
"ES12"  6    3030855 2381070.3
"ES13"  6    3185051 2338763.5
"ES21"  2    3295558   2301408
"ES22"  2    3367609 2247390.8
"ES23"  0    3290383   2215838
"ES24"  1    3430578 2108780.3
"ES30"  2    3159730   2038493
"ES41"  0    3095706   2192837
"ES42"  0    3203476 1927800.5
"ES43"  0    2929076 1937215.6
"ES51"  0    3615668 2116901.3
"ES52"  1    3409599 1874935.8
"ES53"  0    3709127   1860123
"ES61"  0    3030216   1722690
"ES62"  0  3309194.5 1734055.3
"ES70"  0  1800390.3 1016261.3
"FI1C1" 1    5007046   4222733
"FI1D2" 2    5195955   4559601
"FI1D5" 1    5017250   4570523
"FI1D9" 1    5100395 4740287.5
"FR1"   6  3770310.5 2872009.5
"FRB"   3    3695405   2742625
"FRC"   1    3927996   2693985
"FRD"   2    3601146 2938269.5
"FRE"   6    3803574   3009204
"FRF"   4  3998597.5   2851375
"FRG"   2    3507977 2766441.5
"FRH"   0    3371326 2868089.5
"FRI"   2    3552162   2503978
"FRJ"   0    3687057   2321472
"FRK"   2    3894650   2505264
"FRL"   1  4003699.5 2324502.5
"FRM"   0  4246797.5 2116740.3
"FRY1"  0   -2604120   2703846
"FRY2"  6   -2669187   2507016
"FRY3"  0 -2638118.5  779264.4
"FRY4"  0    9991020  -3050530
"FRY5"  0    8728395  -2780427
"HR"    1    4825509 2457686.3
"HU"    3    5031611 2718008.5
"IE"    2  3121504.5   3490938
"IT"    2  4502339.5 2181192.8
"ITC1"  3    4156950 2440637.3
"ITC2"  2    4117523   2516596
"ITC3"  1    4217098 2351550.3
"ITC4"  2    4302959   2500767
"ITF1"  2    4640308 2132906.8
"ITF2"  0    4704800 2076411.5
"ITF3"  0    4730917   1986481
"ITF4"  1  4880878.5 2011984.6
"ITF5"  0  4838518.5   1954352
"ITF6"  0    4872946   1798236
"ITG1"  1  4690059.5   1622663
"ITG2"  0  4237801.5   1888152
"ITH1"  0  4429322.5   2621478
"ITH2"  0    4407685 2558676.5
"ITH3"  0  4464466.5   2507284
"ITH4"  2    4556965   2565468
"ITH5"  2  4404050.5 2380241.5
"ITI1"  6    4412582 2260659.3
end

Thank you for your help!

[Andrew Musau]

spmap is from SSC, as you are asked to explain (FAQ Advice #12). You will want to restrict the displayed coordinates below, but here is some technique.

copy https://gisco-services.ec.europa.eu/distribution/v2/nuts/shp/NUTS_RG_03M_2016_3035_LEVL_2.shp.zip NUTS_RG_03M_2016_3035_LEVL_2.shp.zip, replace
unzipfile NUTS_RG_03M_2016_3035_LEVL_2.shp.zip, replace
spshape2dta NUTS_RG_03M_2016_3035_LEVL_2.shp, saving(euro) replace

* Example generated by -dataex-. For more info, type help dataex
clear
input str5 NUTS_ID float(total x_cc y_cc)
"AT"    2  4631152.5   2728758
"BE1"   1    3924979   3095776
"BE3"   2    3965977   3034637
"BG"    6    5561655 2310162.5
"CY"    0    6426424   1653395
"CZ"    6  4705460.5   2973015
"DE1"   2    4250334   2825138
"DE2"   2  4425542.5   2870462
"DE3"   2  4551963.5 3271179.5
"DE4"   2  4552249.5   3267538
"DE5"   0    4236640   3345796
"DE6"   2  4321787.5   3382180
"DE7"   2  4252034.5   3054615
"DE8"   2    4488061 3408451.5
"DE9"   2    4264924   3295466
"DEA"   0    4151113   3154747
"DEB"   1  4138005.5   2980828
"DEC"   0  4099770.5   2923665
"DED"   2  4555406.5   3109965
"DEE"   2    4437922   3212359
"DEF"   1    4309255   3452599
"DEG"   1    4393193   3088532
"DK"    2  4324183.5   3650497
"EE"    2    5216204   4050577
"EL"    0    5446546   1866405
"EL30"  6    5519687 1747141.6
"EL41"  0    5718016   1894183
"EL42"  0    5781137 1675096.3
"EL43"  1    5677244   1484150
"EL51"  0    5587258 2134529.8
"EL52"  2    5414566   2054253
"EL53"  0    5297153 1992473.8
"EL54"  1    5248981 1900077.8
"EL61"  0    5372253   1911631
"EL62"  0    5232647   1788519
"EL63"  0    5336454   1765899
"EL64"  0    5439385 1824763.5
"EL65"  0    5428266 1675550.5
"ES"    2    3155919   2006099
"ES11"  0  2866097.5 2357738.8
"ES12"  6    3030855 2381070.3
"ES13"  6    3185051 2338763.5
"ES21"  2    3295558   2301408
"ES22"  2    3367609 2247390.8
"ES23"  0    3290383   2215838
"ES24"  1    3430578 2108780.3
"ES30"  2    3159730   2038493
"ES41"  0    3095706   2192837
"ES42"  0    3203476 1927800.5
"ES43"  0    2929076 1937215.6
"ES51"  0    3615668 2116901.3
"ES52"  1    3409599 1874935.8
"ES53"  0    3709127   1860123
"ES61"  0    3030216   1722690
"ES62"  0  3309194.5 1734055.3
"ES70"  0  1800390.3 1016261.3
"FI1C1" 1    5007046   4222733
"FI1D2" 2    5195955   4559601
"FI1D5" 1    5017250   4570523
"FI1D9" 1    5100395 4740287.5
"FR1"   6  3770310.5 2872009.5
"FRB"   3    3695405   2742625
"FRC"   1    3927996   2693985
"FRD"   2    3601146 2938269.5
"FRE"   6    3803574   3009204
"FRF"   4  3998597.5   2851375
"FRG"   2    3507977 2766441.5
"FRH"   0    3371326 2868089.5
"FRI"   2    3552162   2503978
"FRJ"   0    3687057   2321472
"FRK"   2    3894650   2505264
"FRL"   1  4003699.5 2324502.5
"FRM"   0  4246797.5 2116740.3
"FRY1"  0   -2604120   2703846
"FRY2"  6   -2669187   2507016
"FRY3"  0 -2638118.5  779264.4
"FRY4"  0    9991020  -3050530
"FRY5"  0    8728395  -2780427
"HR"    1    4825509 2457686.3
"HU"    3    5031611 2718008.5
"IE"    2  3121504.5   3490938
"IT"    2  4502339.5 2181192.8
"ITC1"  3    4156950 2440637.3
"ITC2"  2    4117523   2516596
"ITC3"  1    4217098 2351550.3
"ITC4"  2    4302959   2500767
"ITF1"  2    4640308 2132906.8
"ITF2"  0    4704800 2076411.5
"ITF3"  0    4730917   1986481
"ITF4"  1  4880878.5 2011984.6
"ITF5"  0  4838518.5   1954352
"ITF6"  0    4872946   1798236
"ITG1"  1  4690059.5   1622663
"ITG2"  0  4237801.5   1888152
"ITH1"  0  4429322.5   2621478
"ITH2"  0    4407685 2558676.5
"ITH3"  0  4464466.5   2507284
"ITH4"  2    4556965   2565468
"ITH5"  2  4404050.5 2380241.5
"ITI1"  6    4412582 2260659.3
end

tempfile mydata
save `mydata'

use euro.dta, clear
merge m:1 NUTS_ID using `mydata', keep(match) nogen
spmap using euro_shp.dta, id(_ID) ndocolor(black) fcolor(none) leg(off)

[Ema Davies]

Thank you, Andrew. Here is what I tried, but it still does not show all the values:

ssc install shp2dta
ssc install mergepoly
ssc install spmap

*  1.  Get  shapefile  

foreach val in 0 1 2 3 {
copy https://gisco-services.ec.europa.eu/distribution/v2/nuts/shp/NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip "NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip"
unzipfile "NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip", replace
}

*  2.  Convert  shapefile  to  Stata  attribute  and  coordinate  datasets  

foreach val in 0 1 2 3 {
shp2dta using  "NUTS_RG_03M_2016_3035_LEVL_`val'",  data("reg--attr_`val'.dta")  coord("reg--coord_`val'.dta")      ///
      genid(stid)  gencentroids(cc) replace
}

*Append the files

use "reg--attr_0.dta", clear
gen NUTSlevel=0
foreach val in 1 2 3 {
append using "reg--attr_`val'.dta"
replace NUTSlevel=`val' if NUTSlevel==.
}
drop stid
gen _ID=_n
order _ID
save "reg--attr.dta", replace

*The data:

* Example generated by -dataex-. For more info, type help dataex
clear
input str5 NUTS_ID float total
"AT"    2
"BE1"   1
"BE3"   2
"BG"    6
"CY"    0
"CZ"    6
"DE1"   2
"DE2"   2
"DE3"   2
"DE4"   2
"DE5"   0
"DE6"   2
"DE7"   2
"DE8"   2
"DE9"   2
"DEA"   0
"DEB"   1
"DEC"   0
"DED"   2
"DEE"   2
"DEF"   1
"DEG"   1
"DK"    2
"EE"    2
"EL"    0
"EL30"  6
"EL41"  0
"EL42"  0
"EL43"  1
"EL51"  0
"EL52"  2
"EL53"  0
"EL54"  1
"EL61"  0
"EL62"  0
"EL63"  0
"EL64"  0
"EL65"  0
"ES"    2
"ES11"  0
"ES12"  6
"ES13"  6
"ES21"  2
"ES22"  2
"ES23"  0
"ES24"  1
"ES30"  2
"ES41"  0
"ES42"  0
"ES43"  0
"ES51"  0
"ES52"  1
"ES53"  0
"ES61"  0
"ES62"  0
"ES70"  0
"FI193" 0
"FI194" 0
"FI195" 0
"FI196" 0
"FI197" 0
"FI1B1" 0
"FI1C1" 1
"FI1C2" 0
"FI1C3" 0
"FI1C4" 0
"FI1C5" 0
"FI1D1" 0
"FI1D2" 2
"FI1D3" 0
"FI1D5" 1
"FI1D7" 0
"FI1D8" 0
"FI1D9" 1
"FR1"   6
"FRB"   3
"FRC"   1
"FRD"   2
"FRE"   6
"FRF"   4
"FRG"   2
"FRH"   0
"FRI"   2
"FRJ"   0
"FRK"   2
"FRL"   1
"FRM"   0
"FRY1"  0
"FRY2"  6
"FRY3"  0
"FRY4"  0
"FRY5"  0
"HR"    1
"HU"    3
"IE"    2
"IT"    2
"ITC1"  3
"ITC2"  2
"ITC3"  1
"ITC4"  2
end

save mydata, replace

*The map:

use mydata, clear

merge 1:1 NUTS_ID using "reg--attr.dta"
drop if _m==2
drop _m
    
spmap total using "reg--coord_0.dta",  id(_ID)  fcolor(Blues) ///
    clmethod(custom) clbreaks(0 1 2 3 4 5 6) ///
    ocolor(white) osize(thin) ///
    legend(label(1 "No") label(2 "1") label(3 "2") label(4 "3") label(5 "4") label(6 "5") label(7 "High")) ///
    polygon (data("reg--coord_2.dta") ocolor(purple) osize(thin))

I would appreciate any other suggestions you might have.

[Andrew Musau]

This is what you get if you tabulate the variable total:

. tab total

      total |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |         45       45.00       45.00
          1 |         15       15.00       60.00
          2 |         28       28.00       88.00
          3 |          3        3.00       91.00
          4 |          1        1.00       92.00
          6 |          8        8.00      100.00
------------+-----------------------------------
      Total |        100      100.00

There are two issues. First, because not all values between 0 and 6 are observed, an ordinal scale will differ from a continuous scale. I will assume a continuous scale.

Second, since you are displaying the map in its entirety, there is a problem when it comes to distinguishing between areas with zero totals within your sample and areas that are outside your sample (i.e., zero observations). If this distinction is important, you need to apply shading to the in-sample zero-observation areas that differs from the default no-shading used for areas outside your sample.

Regarding the construction of the scale, you can define a separate variable for this purpose. This is illustrated below.

foreach val in 0 1 2 3 {
copy https://gisco-services.ec.europa.eu/distribution/v2/nuts/shp/NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip "NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip", replace
unzipfile "NUTS_RG_03M_2016_3035_LEVL_`val'.shp.zip", replace
}

*  2.  Convert  shapefile  to  Stata  attribute  and  coordinate  datasets  

foreach val in 0 1 2 3 {
shp2dta using  "NUTS_RG_03M_2016_3035_LEVL_`val'",  data("reg--attr_`val'.dta")  coord("reg--coord_`val'.dta")      ///
      genid(stid)  gencentroids(cc) replace
}

* Example generated by -dataex-. For more info, type help dataex
clear
input str5 NUTS_ID float total
"AT"    2
"BE1"   1
"BE3"   2
"BG"    6
"CY"    0
"CZ"    6
"DE1"   2
"DE2"   2
"DE3"   2
"DE4"   2
"DE5"   0
"DE6"   2
"DE7"   2
"DE8"   2
"DE9"   2
"DEA"   0
"DEB"   1
"DEC"   0
"DED"   2
"DEE"   2
"DEF"   1
"DEG"   1
"DK"    2
"EE"    2
"EL"    0
"EL30"  6
"EL41"  0
"EL42"  0
"EL43"  1
"EL51"  0
"EL52"  2
"EL53"  0
"EL54"  1
"EL61"  0
"EL62"  0
"EL63"  0
"EL64"  0
"EL65"  0
"ES"    2
"ES11"  0
"ES12"  6
"ES13"  6
"ES21"  2
"ES22"  2
"ES23"  0
"ES24"  1
"ES30"  2
"ES41"  0
"ES42"  0
"ES43"  0
"ES51"  0
"ES52"  1
"ES53"  0
"ES61"  0
"ES62"  0
"ES70"  0
"FI193" 0
"FI194" 0
"FI195" 0
"FI196" 0
"FI197" 0
"FI1B1" 0
"FI1C1" 1
"FI1C2" 0
"FI1C3" 0
"FI1C4" 0
"FI1C5" 0
"FI1D1" 0
"FI1D2" 2
"FI1D3" 0
"FI1D5" 1
"FI1D7" 0
"FI1D8" 0
"FI1D9" 1
"FR1"   6
"FRB"   3
"FRC"   1
"FRD"   2
"FRE"   6
"FRF"   4
"FRG"   2
"FRH"   0
"FRI"   2
"FRJ"   0
"FRK"   2
"FRL"   1
"FRM"   0
"FRY1"  0
"FRY2"  6
"FRY3"  0
"FRY4"  0
"FRY5"  0
"HR"    1
"HU"    3
"IE"    2
"IT"    2
"ITC1"  3
"ITC2"  2
"ITC3"  1
"ITC4"  2
end

gen tot= total+2
save mydata, replace

*Append the files

use "reg--attr_0.dta", clear
gen NUTSlevel=0
foreach val in 1 2 3 {
append using "reg--attr_`val'.dta"
replace NUTSlevel=`val' if NUTSlevel==.
}
drop stid
gen _ID=_n
order _ID
save "reg--attr.dta", replace

use mydata, clear

merge 1:1 NUTS_ID using "reg--attr.dta"
drop if _m==2
drop _m
    
spmap tot using "reg--coord_0.dta",  id(_ID) ///
 fcolor(red*.2 blue*.2 blue*.3 blue*.5 blue*1.0 blue*1.5 blue*2.0 ) ///
    clmethod(custom) clbreaks(1 2 3 4 5 6 7 8) ///
    ocolor(white) osize(thin) ///
    leg(lab(2 "No") lab(3 "1") lab(4 "2") lab(5 "3") lab(6 "4") lab(7 "5") lab(8 "High")) ///
    polygon (data("reg--coord_2.dta") ocolor(purple) osize(thin))

[Ema Davies]

Thank you very much, Andrew, for your help! My concern is that, for some countries present in the sample, such as France, where total > 0, it still shows on the map as if there is no data. If I try using "reg---coord_1.dta" instead of "reg---coord_0.dta" in the code below, it partly resolves the issue, but the problem still persists for some other countries.

 spmap tot using "reg--coord_1.dta",  id(_ID) ///
 fcolor(red*.2 blue*.2 blue*.3 blue*.5 blue*1.0 blue*1.5 blue*2.0 ) ///  
 clmethod(custom) clbreaks(1 2 3 4 5 6 7 8) ///  
 ocolor(white) osize(thin) ///  
 leg(lab(2 "No") lab(3 "1") lab(4 "2") lab(5 "3") lab(6 "4") lab(7 "5") lab(8 "High")) ///
 polygon (data("reg--coord_2.dta") ocolor(purple) osize(thin))

[Andrew Musau]

The issue you are encountering is likely due to a mismatch between the NUTS level in your data and the NUTS level in the coordinate file used to draw the map. Your dataset contains regions at different levels of aggregation (for example, some country-level codes like "FR" and some regional codes such as "FRJ" or "FRY1"). If you then merge this data with a coordinate file that is defined at a specific NUTS level (say NUTS 1 or NUTS 2), any observation that does not exactly match that level will fail to merge. When this happens, the map polygon still exists, but your variable (e.g., total) is missing for that unit, so it is displayed as if there were no data. Switching from coord_0 to coord_1 partially resolves the problem because you are changing the aggregation level of the map layer, but if your data still contain mixed NUTS levels, some mismatches will remain. The solution is to ensure that your data and the coordinate file are defined at the same NUTS level throughout, either by aggregating your data to a single level or by using the coordinate file that exactly corresponds to the level of your NUTS identifiers. The classification is defined here: https://ec.europa.eu/eurostat/web/nuts. However, a quick rule of thumb is to check the length of the code, where:

Length = 2 → NUTS 0
Length = 3 → NUTS 1
Length = 4 → NUTS 2
Length = 5 → NUTS 3