将比例尺（以米为单位）添加到 Stata 中 GPS 坐标的散点图
Add scale bar (in meters) to scatterplot of GPS coordinates in Stata
我正在 Stata 中使用 twoway scatter 绘制一系列 GPS 点。这只是为了可视化这些点的位置以进行一些快速验证。我想为 x 轴和 y 轴添加一个比例尺，以了解这些点实际相距多远。
下面是一些代码，用于输入一些坐标并绘制散点图，指定轴范围：
* Input example data
clear
input key lat lon
1 -17.01639 36.11884
2 -17.01669 36.11916
3 -17.01692 36.11666
4 -17.01575 36.11607
5 -17.01551 36.11619
6 -17.01543 36.11665
7 -17.01581 36.11706
8 -17.01665 36.11672
9 -17.01543 36.11612
10 -17.01724 36.11668
end
* Install Geodist (we will need this later)
cap which geodist
if _rc == 111 ssc install geodist // Installs geodist if necessary
* Create y/x axis ranges using min and max latitude and longitude
sort lat
local ymin = lat[1]
local ymax = lat[_N]
local ydelta = (`ymax' - `ymin')/3
sort lon
local xmin = lon[1]
local xmax = lon[_N]
local xdelta = (`xmax' - `xmin')/3
* Graph points
twoway scatter lat lon, ///
ylabel(`ymin'(`ydelta')`ymax', angle(forty_five)) ///
xlabel(`xmin'(`xdelta')`xmax')
你应该得到一个可爱的基本散点图。这很好，但我不知道这些点之间到底有多远。现在，我们可以使用以下方法了解最小值和最大值 latitude/longitude 之间的距离：
* Store x and y ranges in locals
geodist `ymin' `xmin' `ymax' `xmin'
local yrange = `r(distance)'
geodist `ymin' `xmin' `ymin' `xmax'
local xrange = `r(distance)'
对于此数据，y 跨度为 0.2km，x 跨度为 0.329km。
注意：我知道在地球的球面投影上，ymin 和 ymax 之间的距离可能无法完美地映射到该距离的子集（例如，从 ymin 到 ymax/2 的距离可能不等于 ymin 和 ymax 之间距离的一半）。但是对于我需要这种可视化的目的，以及这些数据通常聚集在彼此相距 2 公里以内的事实，它足够准确。
那么，我如何为每个 x 和 y 创建一个比例尺来显示，例如 100m 的样子？在此示例中，y 比例看起来与 x 比例不同，所以也许我需要找到一种方法来更改 y 范围和 x 范围，以便它们可以表示相同的比例？
我不确定这是否正是你想要的，我在 Stata 中没有任何使用 geo2xy 的经验，但是使用常规 twoway 和一些代数你可以得到近似值这些尺度的值（特别是因为你说你在一个小尺度上工作，近似是可以的）
我把代码放在下面，但是有点长，所以这里是执行摘要：
首先，我们要计算 1 米距离的纬度和经度值。我们可以这样做，因为我们有以公里为单位的范围以及纬度和经度的最小值和最大值。
然后我将这个距离放入一个 twoway connect 图表中，放在我们之前的散点图上。
要获取 1 米距离的值，我执行以下操作：
我们可以假设这几乎是线性的（0 映射到 0 公里，xmax-xmin 映射到 xrange 公里）所以为了得到一米，我们可以将其代入线性方程 y=mx+b，注意b是0，我们的slope = xrange/(xmax-xmin)。然后，我们可以设置1 meter = 1/1000 * 1/slope或1 meter = (xmax-xmin)/(xrange*1000)。在我的示例中，我使用 100 米，因为 1 米太小了，但是要更改比例，您只需将 1/1000 更改为您想要的任何比例。
然后我们可以对 y 轴做同样的事情并将它们添加到数据集并绘制它们。
这是我在 运行 你上面的代码之后用来玩弄它的代码：
local xmeter = (`xmax'-`xmin')/(`xrange'*10)
local ymeter = (`ymax'-`ymin')/(`yrange'*10)
local xbit = (`xmax'-`xmin')/20
local ybit = (`ymax'-`ymin')/20
local xmin_and_bit = `xmin' + `xbit'
local ymin_and_bit = `ymin' + `ybit'
local xmin_and_bit_and_one_meter = `xmin' + `xbit' + `xmeter'
local ymin_and_bit_and_one_meter = `ymin' + `ybit' + `ymeter'
// generate a dummy variable for data vs lengths
gen scale = 0
// add four observations for our distance markers
local new_obs_num = _N + 4
set obs `new_obs_num'
// point for ymin and xmin+bit
replace scale = 1 if _n == _N - 3
replace lat = `ymin' if _n == _N - 3
replace lon = `xmin_and_bit' if _n == _N - 3
// point for ymin and xmin+bit+1
replace scale = 1 if _n == _N - 2
replace lat = `ymin' if _n == _N - 2
replace lon = `xmin_and_bit_and_one_meter' if _n == _N - 2
// point for ymin+bit and xmin
replace scale = 2 if _n == _N - 1
replace lat = `ymin_and_bit' if _n == _N - 1
replace lon = `xmin' if _n == _N - 1
// point for ymin+bit+1 and xmin
replace scale = 2 if _n == _N
replace lat = `ymin_and_bit_and_one_meter' if _n == _N
replace lon = `xmin' if _n == _N
// now graph the scatter plot if scale == 0,
// the x 100 meters if scale == 1,
// and y 100 meters if scale == 2
twoway (scatter lat lon if scale == 0) ///
(connected lat lon if scale == 1, ///
msymbol(+) lcolor(black) mcolor(black)) ///
(connected lat lon if scale == 2, ///
msymbol(+) lcolor(black) mcolor(black)), ///
ylabel(`ymin'(`ydelta')`ymax', angle(forty_five)) ///
xlabel(`xmin'(`xdelta')`xmax') ///
legend(off) ///
note("Points are points and lines show 100 meters")
如果您要创建带有地理坐标的地图，您首先需要将它们转换为平面 (x,y) 坐标。使用 geo2xy（来自 SSC），您可以使用与 Google 地图相同的投影，因此您创建的地图将具有与您在 Google 地图上观察到的完全相同的比例。由于您需要基于距离单位的坐标，因此我选择了真正的墨卡托投影，它以米为单位创建 (x,y) 坐标。默认情况下，投影以中经为中心，这意味着 x 坐标将以距地图中心的米为单位。 y 坐标以距赤道的米为单位。要将它们集中在地图的中间，您只需将它们从平均值中偏移即可。
在 Stata 中创建地图时，您需要确保创建具有正确比例的地图区域（xsize() 和 ysize() 与距离成正比）。
* Input example data
clear
input key lat lon
1 -17.01639 36.11884
2 -17.01669 36.11916
3 -17.01692 36.11666
4 -17.01575 36.11607
5 -17.01551 36.11619
6 -17.01543 36.11665
7 -17.01581 36.11706
8 -17.01665 36.11672
9 -17.01543 36.11612
10 -17.01724 36.11668
end
* project lat/lon to (x,y) in meters using a Mercator projection
geo2xy lat lon, gen(ylat xlon) project(mercator)
* show the projection details and compute the plot's height
return list
local yheight = 6 * `r(aspect)'
* the projection is in meters, centered on mid-longitude, center y coor as well
sum ylat, meanonly
replace ylat = ylat - r(mean)
scatter ylat xlon, ///
xsize(6) ysize(`yheight') ///
ylabel(minmax, nogrid) yscale(off) ///
xlabel(minmax, nogrid) xscale(off) ///
plotregion(margin(small)) graphregion(margin(small)) ///
legend(off) name(mercator_us, replace)
graph export my_map.png, width(600) replace
这是生成的地图：
如果要显示轴单位（以米为单位），您可以删除 xlabel() 和 ylabel() 行。如果你这样做，Stata 可能会重新缩放每个轴以适应合适的间隔，因此会稍微改变比例。如果我让 Stata 决定标签，地图会是这样的：
我正在 Stata 中使用 twoway scatter 绘制一系列 GPS 点。这只是为了可视化这些点的位置以进行一些快速验证。我想为 x 轴和 y 轴添加一个比例尺，以了解这些点实际相距多远。
下面是一些代码，用于输入一些坐标并绘制散点图，指定轴范围：
* Input example data
clear
input key lat lon
1 -17.01639 36.11884
2 -17.01669 36.11916
3 -17.01692 36.11666
4 -17.01575 36.11607
5 -17.01551 36.11619
6 -17.01543 36.11665
7 -17.01581 36.11706
8 -17.01665 36.11672
9 -17.01543 36.11612
10 -17.01724 36.11668
end
* Install Geodist (we will need this later)
cap which geodist
if _rc == 111 ssc install geodist // Installs geodist if necessary
* Create y/x axis ranges using min and max latitude and longitude
sort lat
local ymin = lat[1]
local ymax = lat[_N]
local ydelta = (`ymax' - `ymin')/3
sort lon
local xmin = lon[1]
local xmax = lon[_N]
local xdelta = (`xmax' - `xmin')/3
* Graph points
twoway scatter lat lon, ///
ylabel(`ymin'(`ydelta')`ymax', angle(forty_five)) ///
xlabel(`xmin'(`xdelta')`xmax')
你应该得到一个可爱的基本散点图。这很好，但我不知道这些点之间到底有多远。现在，我们可以使用以下方法了解最小值和最大值 latitude/longitude 之间的距离：
* Store x and y ranges in locals
geodist `ymin' `xmin' `ymax' `xmin'
local yrange = `r(distance)'
geodist `ymin' `xmin' `ymin' `xmax'
local xrange = `r(distance)'
对于此数据，y 跨度为 0.2km，x 跨度为 0.329km。
注意：我知道在地球的球面投影上，ymin 和 ymax 之间的距离可能无法完美地映射到该距离的子集（例如，从 ymin 到 ymax/2 的距离可能不等于 ymin 和 ymax 之间距离的一半）。但是对于我需要这种可视化的目的，以及这些数据通常聚集在彼此相距 2 公里以内的事实，它足够准确。
那么，我如何为每个 x 和 y 创建一个比例尺来显示，例如 100m 的样子？在此示例中，y 比例看起来与 x 比例不同，所以也许我需要找到一种方法来更改 y 范围和 x 范围，以便它们可以表示相同的比例？
我不确定这是否正是你想要的，我在 Stata 中没有任何使用 geo2xy 的经验，但是使用常规 twoway 和一些代数你可以得到近似值这些尺度的值（特别是因为你说你在一个小尺度上工作，近似是可以的）
我把代码放在下面，但是有点长，所以这里是执行摘要：
首先，我们要计算 1 米距离的纬度和经度值。我们可以这样做，因为我们有以公里为单位的范围以及纬度和经度的最小值和最大值。
然后我将这个距离放入一个 twoway connect 图表中，放在我们之前的散点图上。
要获取 1 米距离的值，我执行以下操作：
我们可以假设这几乎是线性的（0 映射到 0 公里，xmax-xmin 映射到 xrange 公里）所以为了得到一米，我们可以将其代入线性方程 y=mx+b，注意b是0，我们的slope = xrange/(xmax-xmin)。然后，我们可以设置1 meter = 1/1000 * 1/slope或1 meter = (xmax-xmin)/(xrange*1000)。在我的示例中，我使用 100 米，因为 1 米太小了，但是要更改比例，您只需将 1/1000 更改为您想要的任何比例。
然后我们可以对 y 轴做同样的事情并将它们添加到数据集并绘制它们。
这是我在 运行 你上面的代码之后用来玩弄它的代码：
local xmeter = (`xmax'-`xmin')/(`xrange'*10)
local ymeter = (`ymax'-`ymin')/(`yrange'*10)
local xbit = (`xmax'-`xmin')/20
local ybit = (`ymax'-`ymin')/20
local xmin_and_bit = `xmin' + `xbit'
local ymin_and_bit = `ymin' + `ybit'
local xmin_and_bit_and_one_meter = `xmin' + `xbit' + `xmeter'
local ymin_and_bit_and_one_meter = `ymin' + `ybit' + `ymeter'
// generate a dummy variable for data vs lengths
gen scale = 0
// add four observations for our distance markers
local new_obs_num = _N + 4
set obs `new_obs_num'
// point for ymin and xmin+bit
replace scale = 1 if _n == _N - 3
replace lat = `ymin' if _n == _N - 3
replace lon = `xmin_and_bit' if _n == _N - 3
// point for ymin and xmin+bit+1
replace scale = 1 if _n == _N - 2
replace lat = `ymin' if _n == _N - 2
replace lon = `xmin_and_bit_and_one_meter' if _n == _N - 2
// point for ymin+bit and xmin
replace scale = 2 if _n == _N - 1
replace lat = `ymin_and_bit' if _n == _N - 1
replace lon = `xmin' if _n == _N - 1
// point for ymin+bit+1 and xmin
replace scale = 2 if _n == _N
replace lat = `ymin_and_bit_and_one_meter' if _n == _N
replace lon = `xmin' if _n == _N
// now graph the scatter plot if scale == 0,
// the x 100 meters if scale == 1,
// and y 100 meters if scale == 2
twoway (scatter lat lon if scale == 0) ///
(connected lat lon if scale == 1, ///
msymbol(+) lcolor(black) mcolor(black)) ///
(connected lat lon if scale == 2, ///
msymbol(+) lcolor(black) mcolor(black)), ///
ylabel(`ymin'(`ydelta')`ymax', angle(forty_five)) ///
xlabel(`xmin'(`xdelta')`xmax') ///
legend(off) ///
note("Points are points and lines show 100 meters")
如果您要创建带有地理坐标的地图，您首先需要将它们转换为平面 (x,y) 坐标。使用 geo2xy（来自 SSC），您可以使用与 Google 地图相同的投影，因此您创建的地图将具有与您在 Google 地图上观察到的完全相同的比例。由于您需要基于距离单位的坐标，因此我选择了真正的墨卡托投影，它以米为单位创建 (x,y) 坐标。默认情况下，投影以中经为中心，这意味着 x 坐标将以距地图中心的米为单位。 y 坐标以距赤道的米为单位。要将它们集中在地图的中间，您只需将它们从平均值中偏移即可。
在 Stata 中创建地图时，您需要确保创建具有正确比例的地图区域（xsize() 和 ysize() 与距离成正比）。
* Input example data
clear
input key lat lon
1 -17.01639 36.11884
2 -17.01669 36.11916
3 -17.01692 36.11666
4 -17.01575 36.11607
5 -17.01551 36.11619
6 -17.01543 36.11665
7 -17.01581 36.11706
8 -17.01665 36.11672
9 -17.01543 36.11612
10 -17.01724 36.11668
end
* project lat/lon to (x,y) in meters using a Mercator projection
geo2xy lat lon, gen(ylat xlon) project(mercator)
* show the projection details and compute the plot's height
return list
local yheight = 6 * `r(aspect)'
* the projection is in meters, centered on mid-longitude, center y coor as well
sum ylat, meanonly
replace ylat = ylat - r(mean)
scatter ylat xlon, ///
xsize(6) ysize(`yheight') ///
ylabel(minmax, nogrid) yscale(off) ///
xlabel(minmax, nogrid) xscale(off) ///
plotregion(margin(small)) graphregion(margin(small)) ///
legend(off) name(mercator_us, replace)
graph export my_map.png, width(600) replace
这是生成的地图：
如果要显示轴单位（以米为单位），您可以删除 xlabel() 和 ylabel() 行。如果你这样做，Stata 可能会重新缩放每个轴以适应合适的间隔，因此会稍微改变比例。如果我让 Stata 决定标签，地图会是这样的：