标签: shapefile

我有一个ESRI shapefile(来自这里:http://pubs.usgs.gov/ds/425/).我希望使用python以给定的纬度/经度从形状文件(在这种情况下为表面材料)中查找信息.

解决这个问题的最佳方法是什么？

谢谢.

最终解决方案

#!/usr/bin/python

from osgeo import ogr, osr

dataset = ogr.Open('./USGS_DS_425_SHAPES/Surficial_materials.shp')
layer = dataset.GetLayerByIndex(0)
layer.ResetReading()

# Location for New Orleans: 29.98 N, -90.25 E
point = ogr.CreateGeometryFromWkt("POINT(-90.25 29.98)")

# Transform the point into the specified coordinate system from WGS84
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)
coordTransform = osr.CoordinateTransformation(
        spatialRef, layer.GetSpatialRef())

point.Transform(coordTransform)

for feature in layer:
    if feature.GetGeometryRef().Contains(point):
        break

for i in range(feature.GetFieldCount()):
    print feature.GetField(i)

有没有办法在 R 中打开 DWG 文件并将其另存为 ESRI shapefile？我还不知道有什么可用的软件包。

提前致谢。

我有一个托管在 Windows 2008 Server RT 虚拟机上的 PostgreSQL 数据库（是的，我知道它应该托管在 Linux 虚拟机上，但这是我的组织规定的。唉……）

我们的 GIS 人员将大量 shapefile 转储到存储库中。我们希望有一个自动进程作为计划任务遍历文件夹。我们希望将这些添加到我们目前正在开发的其他一些流程的 Postgres 地理数据库中

我希望浏览大量的形状文件，并将它们的几何形状和文件名加载到数据库中。

这是我迄今为止所工作的摄取功能的核心部分的要点

import os, subprocess
base_dir = r"c:\shape_file_repository"
full_dir = os.walk(base_dir)
shapefile_list = []
for source, dirs, files in full_dir:
    for file_ in files:
        if file_[-3:] == 'shp':
            #print "Found Shapefile"
            shapefile_path = base_dir + '/' + file_
            shapefile_list.append(shapefile_path)
for paths in shapefile_list:  
    #This is the part where I keep running into trouble. os.system also didnt work
    temp_bat = open(r"c:\temp\temp_shp.bat", "w")
    temp_bat.write(r'start /D …

我有莫桑比克铁路的 Shapefile，并使用下面的代码生成了铁路沿线的 100 个随机点。

我的问题很简单，但我找不到答案：如何计算铁路沿线各点之间的距离？

我不需要欧几里得距离，我想要沿着铁轨从 A 点到 B 点的距离。

提前致谢！

library(sp)
library(rgdal)
library(spgrass6)
library(maptools)
library(igraph)
library(fields)

railroads <- readShapeLines("MOZ_rails.shp")

#Generate 100 random points, and put them on a matrix:
RandomPoints<-spsample(railroads, 100, type="random")

我有一个折线形状文件，它代表城市道路网的一部分。我的 shapefile 包含多个线路/街道段（在我的例子中为 58）。通过使用 R-cran，我想进一步将折线段划分为具有相等长度（例如 10 m）的较小部分。

提供一个想法：

当我将折线形状文件导入 R 并创建一个数据框时，它看起来像：

# Import the polyline shapefile into a SpatialLinesDataFrame object:
# library(sp)
sp.roads <- readOGR(dsn="/Users/mgv/Documents/project",layer="roads_sr_corr") 

# Create a dataframe from the SpatialLinesDataFrame 
# library(stplanr)    
linedf <- line2df(sp.roads)
str(linedf)
> str(linedf)
'data.frame':   58 obs. of  4 variables:
 $ fx: num  13.39991 13.40138 13.40606 13.40232 13.40177 ...
 $ fy: num  42.35066 42.35412 42.35599 42.34514 42.34534 ...
 $ tx: num  13.40150 13.40119 13.40591 13.40246 13.40182 ...
 $ ty: num  42.35026 42.35386 42.35602 42.34530 …

使用这里提供的shapefile ,我正在尝试两个合并苏丹和南苏丹的多边形,以便我在2010年获得苏丹的边界.我在R中提供shapefile的代码是

library(rgdal)
africa    <- readOGR(dsn =  "Data/Shapes", layer = "AfricanCountries")
class(africa)

[1] "SpatialPolygonsDataFrame" attr(,"package") [1] "sp"

Run Code Online (Sandbox Code Playgroud)

我试过不同的包和解决方案raster::intersect,rgeos::gIntersect或者maptools::unionSpatialPolygons.我总是得到一个空间对象,它丢失了属于两个多边形或根本没有联合的数据.

到目前为止,我有以下解决方案似乎不是很方便:

# Seperate Polygons and Data
sp_africa  <- SpatialPolygons(africa@polygons, proj4string = CRS(proj4string(africa)))
dat_africa <- africa@data

# Delete South Sudan from dataframe
dat_africa <- dat_africa[-which(dat_africa$COUNTRY == "South Sudan"),]

# Add row.names to data according to polygon ID'S
rownames(dat_africa) <- dat_africa$OBJECTID

# Get all ID's
allIDs <- africa$OBJECTID

# Get the ID for Sudan, Primary …

我希望这不是太琐碎，但我真的找不到答案，而且我对这个话题太陌生，无法自己提出替代方案。所以这里是问题：

我有两个 shapefile x 和 y，它们代表 Sentinel2 卫星图像的不同处理级别。x 包含大约 1.300.000 个多边形/段，完全覆盖了图像扩展，没有任何进一步的重要信息。
y 包含大约 500 个多边形，代表图像的无云区域（也覆盖了图像的大部分，除了一些“云洞”）以及 4 列（传感器、时间...）

我正在尝试将图像信息添加到 x 中 x 被 y 覆盖的位置。很简单？我就是找不到不花几天时间就可以实现的方法。

我将 x 作为一个简单的特征 {sf} 读入，因为使用 shapefile / readOGR 读取它需要很长时间。我用 y 尝试了不同的东西

当我尝试合并（x，y）时，我只能取一个 sf，因为合并不支持两个 sf。合并 x（作为 sf）和 y（作为 shp）给我错误“无法分配大小为 13.0 Gb 的向量”

所以我尝试了sf::st_join(x,y)，它支持两个变量都是 sf 但现在仍然没有完成 28 小时

sf::st_intersect(x,y) 10.000 段子集大约需要 9 分钟，因此对于整个片段来说可能不会快很多。

可以将 x 子集为几个较小的部分来解决整个问题还是有另一个简单的解决方案？我可以对我的工作区做些什么来使合并工作，还是根本没有加入那么多多边形的捷径？

非常感谢，我希望我的描述不会太模糊！

我的小工作站：
win 7 64 位
8 GB RAM
intel i7-4790 @ 3,6 GHz

library(raster)
admin <- getData('GADM', country='FRA', level=2)

set.seed(123)

id <- data.frame(ID_2 = admin@data$ID_2, day1 = sample(1:20,96,replace = T), 
                                   day2 = sample(50:80,96,replace = T), 
                                   day3 = sample(120:140,96,replace = T), 
                                   day4 = sample(200:230,96,replace = T))

 admin.shp <- merge(admin,id)

如果我想用要么色彩的情节day1,day2,day3或者day4,我可以这样做:

plot(admin.shp, col = admin.shp@data$day1)

我想要做的是生成某种从第1天到第230天的.gif文件或动画:

  range <- c(min(id$day1),max(id$day4))

如果日期与值匹配day1,则多边形应变为绿色,如果在第2天匹配多边形应变为蓝色,如果在第3天匹配多边形应变为橙色,如果在第4天匹配多边形应变为红色

如果我必须为一个列(例如day1)执行此操作,我可以这样做:

library(magick)
c(min(id$day1),max(id$day1)) # 1, 20  

for(i in 1:20){

   breaks.pl <- c(0, i, 21) 
   col.pl <- c("green4","white")
   cuts.pl <- cut(data.frame(admin.shp)[, "day1"],breaks = breaks.pl)

   png(paste0(i,".png"), width …

I have a folder with about 100 point shapefiles that are locations obtained while scat sampling of an ungulate species. I would like to merge all these point shapefiles into one shapefile in R. All the point data were in .gpx format initially which I then changed to shapefiles.

我对 R 相当陌生，所以我对如何做到这一点感到非常困惑，并且找不到合并或组合多个 shapefile 的代码。任何建议将不胜感激。谢谢！

例如：[360590, 555610] - [lng, lat] 以米为单位来自谷歌地图 api - GeoJson 数据

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                360590,
                555610
              ],
              [
                360590,
                555555.0128
              ],
              [
                360590,
                555540
              ],
              [
                360592.4439,
                555540
              ],
              [
                360600,
                555540
              ],
              [
                360600,
                555518.8277
              ]
            ]
          ]
        ]
      }
    }
  ]
}

这里，[360590, 555610] - [X, Y] 坐标以米为单位，现在我们需要在谷歌地图上显示这个坐标，有什么解决办法吗？

我们还必须使用 addGeoJson 或 loadGeoJson 方法，因为我们在 GeoJson 文件中有 200MB 数据。现在我们需要在谷歌地图上显示这个坐标，有什么解决办法吗？