Posting hotspots and coldspots points on a map using "spmap" command

[Sumon Khan]

Dear All,
I am working with a DHS data set and calculated hotspots and coldspots of low birth weight in Bangladesh using the Stata "getisord" command. Now my target is to present hotspots and coldspots points on a map so that readers can easily understand- as like a map for a region of the USA regarding sex abuse which is attached herewith. However, when I generated the hotspot and coldspots and going to present them on a map using the "spmap" command, it is not working. I use the following command-

spmap go_z_SBI_yes_b using "administration_1cor", id(id) clm(custom) clb(-100 -2.576 -1.960 1.960 2.576 100) fcolor(ebblue eltblue white orange red) legtitle("{it: z}-value") legstyle(1) legcount legend(size(*1.8))

where go_z_SBI_yes_b is value of hotspots and coldspots inside Bangladesh and its eight regions and "administration_1cor" is the cordinates of Bangladesh. The above command produce a blank map with hotsopot and coldspot value (z-value) value showing at the below of map as a lengend. What I am expecting that this command was worked properly sinec results generated, however, these hotspots and coldspots points are not showing. I appreciate if anybody help me to handle this issue- I want present these hotspots and coldspots points on a map. Thanks in advance!

in advance. Please see the attached.

[Scott Merryman]

I think you need to use the -point()- option within spmap.

[Sumon Khan]

Scott Merryman Many thanks. However, could you please explain a bit more?

[Scott Merryman]

You need to create categories for your z-variable to be used the by() option.

clear
//Shapefile:  https://data.humdata.org/dataset/administrative-boundaries-of-bangladesh-as-of-2015

cd "C:\Users\scott\Desktop\temp\bgd_adm_bbs_20201113_SHP"
spshape2dta bgd_admbnda_adm2_bbs_20201113, replace
spshape2dta bgd_admbnda_adm1_bbs_20201113, replace  

//create random points and variable to plot
use bgd_admbnda_adm2_bbs_20201113
keep  _CX _CY
gen x = _CX + rnormal(0, .01)
gen y = _CY + rnormal(0, .01)
gen z =  floor(5*runiform() + 1)
save temp1, replace

use bgd_admbnda_adm1_bbs_20201113,clear

spmap  using bgd_admbnda_adm1_bbs_20201113_shp, id(_ID) /// 
    point(xcoord(x) ycoord(y)  data(temp1)  by(z) fcolor(Rainbow) legenda(on) )

[Sumon Khan]

Scott Merryman, many thanks for your kind reply- it was really helpful. I have updated the codes and re-run the analysis, however, this produced further difficulty. Map produced is unusual which I attached herewith.

My further issue is you estimated random number ("z" in your calculation) and plotted it in the map. However, in my case z has a specific value (LBW in the attached data file). I need to classify these values as cold spot (if LBW= -100 to -1.96), non-significant (LBW=-1.96 to 1.96), and hot spot (LBW=1.96 to 100) and point these in the map using the "bgd_admbnda_adm1_bbs_20201113" shape as per your previous calculation with three different colors (cold spot, hot spot, non-significant).

Shuld you be kind enough for another look and guide me how can I get this thing done? Many thanks in advance.

[Scott Merryman]

Longitude is the x- coordinate (it has an approximate range from 88 to 93) and latitude is the y-coordinate (with a range 20.5 to 27). See the graph with the x,y grid.

You can use -egen,cut()- to recode the variable:

use test,clear
capture drop spots
egen spots = cut(LBW), at(-100, -1.96, 1.96, 100) icode
label define spot 0 "Cold" 1 "Non-Significant" 2 "Hot", replace
label values spots spot
save,replace

use bgd_admbnda_adm1_bbs_20201113,clear
spmap  using bgd_admbnda_adm1_bbs_20201113_shp, id(_ID) /// 
    point(xcoord(y) ycoord(x)  data(test)  by(spots)  /// 
    fcolor(  blue%50 green%50  red%50) legenda(on) ) /// 
     xlabel(,grid) aspect(1.3) freestyle

[Sumon Khan]

Scott Merryman, many thanks. It is now working. I can finish my research now.

[Sumon Khan]

Scott Merryman , many thanks for your previous reply which helped me a lot in addressing the problem that I had. I am now working with the geographically weighted regression. I have run the model with five explanatory variables and save the regression coefficients in a data file. What I am trying to do now is to plot these coefficients in five different maps (for five different explanatory variables) of Bangladesh, placed side by side as shown in the map attached. Coefficients should be divided into 3 or four equal lengths. I have attached herewith a rough data file that is like the data file I am now working on. I would appreciate your help in addressing this issue.

Many thanks in advance.

[Scott Merryman]

You can use -egen cut- to create the categories.

quietly summarize sur_no

local width = (ceil(r(max)*10)/10 - floor(r(min)*10)/10)/5
display `width'
egen group = cut(sur_no), /// 
    at(`=floor(r(min)*10)/10'(`width')`=ceil(r(max)*10)/10') icodes

save sur_no_data,replace

tabstat sur_no, by(group) stat(min max)
use bgd_admbnda_adm1_bbs_20201113,clear

spmap  using bgd_admbnda_adm1_bbs_20201113_shp, id(_ID) /// 
    point(ycoord(east) xcoord(north)  data(sur_no_data)  by(group)  /// 
    fcolor(Blues) legenda(on) )

[Cae Revoredo]

Hey Scott,

I've used your code, but I have a question related to the use of the by() function. When I want to change the size of the points, only one of them changes. Is there any reason to this and is there any solution? Thanks in advance