让 Rpy2 与 rgdal 一起工作到空间子集点答案

【问题标题】：Getting Rpy2 to work with rgdal to spatially subset points让 Rpy2 与 rgdal 一起工作到空间子集点
【发布时间】：2015-10-10 06:21:55
【问题描述】：

所以我有一些已经可以工作的 R 代码。这会以http://robinlovelace.net/r/2014/07/29/clipping-with-r.html 的方法从数据和空间子集中获取一堆点，以针对 shapefile。

#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。有人知道发生了什么吗？

【问题讨论】：

标签： python r gis rpy2

【解决方案1】：

一种方法是将 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)

【讨论】：