我已经发布了此question on the GIS stack exchange,但是访问量并不大。
我是一个GIS用户,使用R进行统计已经有几年了,很高兴看到发布了许多新的GIS功能(栅格,shapefile,rgdal等)。
我已经在R中完成了许多数据库和表的操作,因此能够从shapefile中添加和删除属性是一种强大的潜力。
我希望我只是想念那里的东西,但是我似乎找不到一种从shapefile属性表添加或删除属性的好方法。
Overflow的任何人都可以回复我的GIS帖子吗?还是可以建议我如何将列表属性表放入数据框并返回以替换当前属性表?
自我的原始帖子以来取得了一些进展:
[这是在弄清楚如何使用属性表appart(.dbf),添加内容,现在我试图将其放回一起以替换原始dbf。
>libary(raster); library(rgdal); library(shapefiles)
>shp<-shapefile(Shape) # D.C. area airport polygons
>summary(shp) #Shapefile properties
Object of class SpatialPointsDataFrame
Coordinates:
min max
coords.x1 281314.2 337904.7
coords.x2 4288867.0 4313507.0
Is projected: TRUE
proj4string :
[+proj=utm +zone=18 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0]
Number of points: 4
Data attributes:
ObjectID LOCID NAME FIELD STATE STATE_FIPS ACAIS TOT_ENP TYPE
Min. :134.0 Length:4 Length:4 Length:4 Length:4 Length:4 Length:4 Min. : 271 Length:4
1st Qu.:242.8 Class :character Class :character Class :character Class :character Class :character Class :character 1st Qu.: 3876 Class :character
Median :339.0 Mode :character Mode :character Mode :character Mode :character Mode :character Mode :character Median :3409113 Mode :character
Mean :310.0 Mean :3717251
3rd Qu.:406.2 3rd Qu.:7122488
Max. :428.0
>shp.AT<-read.dbf(gsub(".shp", ".dbf", Shape), header=TRUE) #Read in the attribute table from the .dbf
> shp.AT # First object in the dbf looks like an attribute table!
$dbf
ObjectID LOCID NAME FIELD STATE STATE_FIPS ACAIS TOT_ENP TYPE
1 134 ADW Andrews AFB <NA> Maryland 24 Y 5078 Military
2 279 DCA Washington National <NA> Virginia 51 Y 6813148 <NA>
3 399 HEF Manassas Regional Harry P Davis Field Virginia 51 Y 271 Regional
4 428 IAD Washington Dulles International <NA> Virginia 51 Y 8050506 International
$header
$header$file.version
[1] 3
$header$file.year
[1] 113
$header$file.month
[1] 4
$header$file.day
[1] 12
$header$num.records
[1] 4
$header$header.length
[1] 321
$header$record.length
[1] 148
$header$fields
NAME TYPE LENGTH DECIMAL
1 ObjectID N 10 0
2 LOCID C 5 0
3 NAME C 45 0
4 FIELD C 30 0
5 STATE C 24 0
6 STATE_FIPS C 2 0
7 ACAIS C 1 0
8 TOT_ENP N 11 0
9 TYPE C 20 0
>shp.tab<-as.data.frame(shp.AT[1]) # Grab the first object of the .dbf as a data.frame
> shp.tab # First list object
dbf.ObjectID dbf.LOCID dbf.NAME dbf.FIELD dbf.STATE dbf.STATE_FIPS dbf.ACAIS dbf.TOT_ENP dbf.TYPE
1 134 ADW Andrews AFB <NA> Maryland 24 Y 5078 Military
2 279 DCA Washington National <NA> Virginia 51 Y 6813148 <NA>
3 399 HEF Manassas Regional Harry P Davis Field Virginia 51 Y 271 Regional
4 428 IAD Washington Dulles International <NA> Virginia 51 Y 8050506 International
> shp.tab$NewAT<-1:nrow(shp.tab) # Add my new attribute
> shp.tab # Added my new attribute, now to get this back into my shapefile
dbf.ObjectID dbf.LOCID dbf.NAME dbf.FIELD dbf.STATE dbf.STATE_FIPS dbf.ACAIS dbf.TOT_ENP dbf.TYPE NewAT
1 134 ADW Andrews AFB <NA> Maryland 24 Y 5078 Military 1
2 279 DCA Washington National <NA> Virginia 51 Y 6813148 <NA> 2
3 399 HEF Manassas Regional Harry P Davis Field Virginia 51 Y 271 Regional 3
4 428 IAD Washington Dulles International <NA> Virginia 51 Y 8050506 International 4
>write.dbf(shp.tab, gsub(".shp", ".dbf", Shape)) # Knew this wouldn't work, but demonstrate attempt to write this object as .dbf.
ERROR:
invalid subscript type 'list'
> shp.AT[1]<-shp.tab # Try replacing the old Object[1] with my new table containing the new attribute.
> shp.AT # The table portion fo the shp.AT is gone. No attributes.
$dbf
[1] 134 279 399 428
$header
$header$file.version
[1] 3
$header$file.year
[1] 113
$header$file.month
[1] 4
$header$file.day
[1] 12
$header$num.records
[1] 4
>write.dbf(shp.AT, gsub(".shp", ".dbf", Shape)) # If I go ahead and overwrite anyway...
我的属性不见了,取而代之的是包含单个字段“ dataframe”的属性表。我递归脚本以再次读取属性表。
> shp.tab
dataframe
1 134
2 279
3 399
4 428
所以,我想我很近。谁能帮我把它重新回到shapefile属性表中?还是有更好的方法?
谢谢
我不确定我是否完全理解您要做什么。看起来您只是想向属性表添加新列?如果这是正确的,则就像对待任何数据框一样对待它。
library(rgdal)
dsn <- system.file("vectors", package = "rgdal")
shp<-readOGR(dsn = dsn, layer = 'cities')
shp$NewAT<-1:nrow(shp)
这与我系统上的shapefile完美配合。我通常使用readOGR()函数依靠rgdal读取shapefile。我相当确定您调用的shapefile()函数也称为rgdal。
已编辑以添加可复制的数据集。
Iam使用R进行GIS监视已经有好几年了,通常与QGIS结合使用。为了处理属性表,我通常将形状另存为CSV,将几何另存为WKT(您可以在QGIS保存对话框中找到这些选项)。在下一步中,我读取R中的所有csv(shapes)并进行统计,连接等。最后,将它们写回到HDD并将其加载回QGIS(无需拖放即可导入对话框)并保存它们作为shapefile。
HTH,乔