python - 让 Rpy2 与 rgdal 一起工作以获取空间子集点

标签 python r gis rpy2

所以我有一些已经可以工作的 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/

相关文章:

python - python中是否有任何采用距离矩阵的好的层次聚类包?

r - 如何合并具有相同列名的多个数据框?

google-maps - 开放 GIS 的标准,如 openstreetmap、cloudmade?

go - 如何将东坐标北坐标转换为经纬度?

r - 使用来自 2 个矩阵的值来索引 R 中的第三个矩阵

python - 如何通过 PyQGIS 使用 cpt-city 目录中的颜色渐变

python - 使用 OpenCV Python 计算相机世界位置

python - 如何将 dict 的值写入 openpyxl 中的空(新)列?

python - 使用 np.einsum 进行光线转换的矢量化范围

r - (R lme4)eval 中的错误(替换(expr),envir,enclos): (maxstephalfit) PIRLS step-halvings failed to reduce deviance in pwrssUpdate