所以我有一些已经可以工作的 R 代码。这采用 http://robinlovelace.net/r/2014/07/29/clipping-with-r.html 的方法从数据和空间子集中针对形状文件获取一堆点。 .
#data is a .csv file with lon lat points
data_points <- SpatialPoints(data)
proj4string(data_points) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
data_ll <- spTransform(data_points, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
melbourne <- readOGR("melbourne_australia.land.coastline","melbourne_australia_land_coast") #this is a shapefile from https://mapzen.com/data/metro-extracts
subset <- data_ll[melbourne,]
情节(墨尔本) 点(子集)
我正在尝试将其转换为相应的 rpy2 脚本。到目前为止我已经;
import pandas as pd
import numpy as np
import rpy2.robjects as ro
import rpy2.robjects.numpy2ri
from rpy2.robjects.packages import importr
rgdal = importr('rgdal')
base = importr('base')
rpy2.robjects.numpy2ri.activate()
data = pd.read_csv('sim.csv')
data = data.values
coordinates = ro.r['coordinates']
proj4string = ro.r['proj4string']
spTransform = ro.r['spTransform']
readOGR = ro.r['readOGR']
SpatialPoints = ro.r['SpatialPoints']
CRS = ro.r['CRS']
class_r = ro.r['class']
key = CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
data_points = SpatialPoints(data, proj4string = key)
data_ll = spTransform(data_points, key)
melbourne = readOGR("melbourne_australia.land.coastline", "melbourne_australia_land_coast")
subset = data_ll[melbourne,]
在最后一行失败,并出现错误 TypeError: 'RS4' object is not subscriptable。有人知道发生了什么事吗?
最佳答案
实现此目的的一种方法是将 R 代码转换为函数,然后将其作为包导入。
这是 R 代码:
library(rgdal)
library(sp)
#import data
data <- read.csv("sim.csv", header = F)
subset_points <- function(data){
data_points <- SpatialPoints(data, proj4string=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
data_ll <- spTransform(data_points, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
sink("/dev/null")
melbourne <- readOGR("melbourne_australia.land.coastline", "melbourne_australia_land_coast")
sink()
subset <- data_ll[melbourne,]
final <- as.data.frame(subset)
return(final)
}
这是 Python 代码:
import rpy2.robjects as ro
import rpy2.robjects.numpy2ri
from rpy2.robjects.packages import importr
from rpy2.robjects.packages import SignatureTranslatedAnonymousPackage
with open('subset_data.R') as fh:
rcode = os.linesep.join(fh.readlines())
subset = SignatureTranslatedAnonymousPackage(rcode, "subset")
rpy2.robjects.numpy2ri.activate()
data = pd.read_csv('sim.csv')
data = data.values
final = subset.subset_points(data)
print(np.array(final).T)
关于python - 让 Rpy2 与 rgdal 一起工作以获取空间子集点,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/33050987/