setwd("/shapefiles/uhf")#shapefile obtained from http://sedac.ciesin.columbia.edu/geonetwork/srv/en/cu-gis_metadata.show?id=3569&currTab=advanced#
`enter code here`map<-readOGR(dsn=".",layer="UHF_42_DOHMH_2009")
`enter code here`plot(map)
`summary(map)` #shows units are in us-ft#
#convert the coordinate system of shapefile to lat/long#
coordinates(map)
proj4string(map)=CRS("+init=esri:102718")
map=spTransform(map,CRS("+proj=longlat +datum=WGS84"))
summary(map)
.
> summary(map) #shows units of measurement are in us-ft Below is the output#
Object of class SpatialPolygonsDataFrame
Coordinates:
min max
x 913090.8 1067310
y 120053.5 272932
Is projected: TRUE
proj4string :
[+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333
+lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
+units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0]
Data attributes:
OBJECTID UHFCODE SHAPE_Leng SHAPE_Area
Min. : 1.0 Min. : 0.0 Min. : 34713 Min. : 35132269
1st Qu.:11.5 1st Qu.:203.5 1st Qu.: 55839 1st Qu.:101843450
Median :22.0 Median :303.0 Median : 73089 Median :164519524
Mean :22.0 Mean :281.7 Mean : 90960 Mean :194568186
3rd Qu.:32.5 3rd Qu.:403.5 3rd Qu.:113901 3rd Qu.:269127674
Max. :43.0 Max. :504.0 Max. :250903 Max. :731130742
UHF_NEIGH BOROUGH
Bayside - Little Neck : 1 Bronx : 7
Bedford Stuyvesant - Crown Heights: 1 Brooklyn :11
Bensonhurst - Bay Ridge : 1 Manhattan :10
Borough Park : 1 N/A : 1
Canarsie - Flatlands : 1 Queens :10
(Other) :37 Staten Island: 4
NA's : 1
#selecting a part of the shapefile#
brooklyn<-subset(map, UHFCODE==205 | UHFCODE==206 | UHFCODE==209,)
str(brooklyn)
summary(brooklyn)
#converting selected shapefile in to owin#
boundary<-as(brooklyn, "owin")
plot(boundary)
summary(boundary) #Below is summary of the output#
Window: polygonal boundary
4 separate polygons (no holes)
vertices area relative.area
polygon 1 16 3.35333e-06 0.000764
polygon 2 205 1.10311e-03 0.251000
polygon 3 173 1.69957e-03 0.387000
polygon 4 178 1.58568e-03 0.361000
enclosing rectangle: [-74.04249, -73.9536] x [40.58028, 40.67145] units
Window area = 0.00439172 square units
#This codes are run after transferring the TB cases (.csv file) into r #
a<-as(brooklyn_tb, "ppp")
str(a)
summary(a)
a$window=boundary plot(a, pch=16)
summary(a$window)
#compute K and L functions#
BROOKLYN<-Kest(a, correction="trans")
plot(BROOKLYN)
我的问题是将 shapefile 转换为自己的对象后显示的“relative.area”是什么？
我必须在这个 shapefile 上绘制 TB 病例的纬度和经度，并在其上运行 Kest() 函数。我需要知道 r 使用什么测量单位来进行 Kest() 计算？有没有办法改变单位？运行 Kest() 函数后，布鲁克林图上的“r”单位是什么？我不能发布图片，因为我没有足够的“声誉”！！！！但这里是潜在 K(r) 图的链接：http ://rgm3.lab.nig.ac.jp/RGM-files/R_CC/result/spatstat/thomas.estK.Rd_001_large.png
