将 shapefile 纬度转换为真实纬度答案

【问题标题】：Converting shapefile latitude to real one将 shapefile 纬度转换为真实纬度
【发布时间】：2017-09-30 12:51:02
【问题描述】：

我有一个地理坐标列表（纬度、经度）和一个具有不同图层的 shapefile。我希望能够确定每个坐标属于哪个图层。

但 shapefile (.shp) 有 poligons，其中纬度和经度以奇怪范围内的数字表示，例如 120724.86008864 和 484497.34058312

我知道 prj 文件包含有关如何进行此转换的信息，但我似乎不明白如何进行。就是这样：

PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199432955] ],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",155000],PARAMETER["False_Northing",463000],PARAMETER["Central_Meridian",5.38763888888889],PARAMETER["Scale_Factor",0.9999079],PARAMETER["Latitude_Of_Origin" ,52.15616055555555],UNIT["米",1]]

具体问题是如何将常规纬度/经度点转换为 shapefile。

在 Python 中工作，使用这个库 http://gdal.org/python/

感谢任何帮助。

【问题讨论】：

您要获取坐标吗？还是解析prj文件？
我可以毫无问题地解析文件，但 shapefile 中的纬度/经度在不同的比例或范围内。例如，纬度 120724.86008864 这显然是错误的。我有正常的纬度，例如 52.3605883。所以我想知道我必须对我的法线坐标应用哪种转换才能像文件的那些。然后我就可以将它们与层中的 poligon 相交。

标签： python shapefile

【解决方案1】：

        import re
        regex = "\d{1,3}\.\d+"
        s ="""PROJCS["RD_New",GEOGCS["GCS_Amersfoort",DATUM["D_Amersfoort",SPHEROID["Bessel_1841",6377397.155,299.1528128]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199432955]],PROJECTION["Double_Stereographic"],PARAMETER["False_Easting",155000],PARAMETER["False_Northing",463000],PARAMETER["Central_Meridian",5.38763888888889],PARAMETER["Scale_Factor",0.9999079],PARAMETER["Latitude_Of_Origin",52.15616055555555],UNIT["Meter",1]] """

        m = re.search(regex, s)

        if m:
            print m.groups()

【讨论】：

你的意图是什么？
我想从字符串中获取小数点后最多 3 位的所有小数。例如：12.55 3.44 321.11

【解决方案2】：

# define input
shape_file = "file.shp"
o_lat = 52.3605883
o_lon = 4.8593157

# geospatial bureocracy
driver = ogr.GetDriverByName('ESRI Shapefile')
shape = driver.Open(shape_file)
layer = shape.GetLayer()
geo_ref = layer.GetSpatialRef()
point_ref = ogr.osr.SpatialReference()
point_ref.ImportFromEPSG(4326)
ctran = ogr.osr.CoordinateTransformation(point_ref, geo_ref)

# critical part: transform longitude/latitude to the shapefile's projection
[t_lon, t_lat, z] = ctran.TransformPoint(o_lon, o_lat)
print('original coords', o_lon, o_lat)
print('transformed coords', t_lon, t_lat)

# create the needle
point = ogr.Geometry(ogr.wkbPoint)
point.SetPoint_2D(0, t_lon, t_lat)
layer.SetSpatialFilter(point)

# look it up
for feature in layer:
    polygon = feature.GetGeometryRef()
    if polygon.Contains(point):
        print('Found it', feature.ExportToJson()

【讨论】：