将 shapefile 纬度转换为真实纬度

Question

我有一个地理坐标列表（纬度、经度）和一个包含不同图层的 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],参数["Scale_Factor",0.9999079],参数["Latitude_Of_Origin",52.15616055555555],单位["Meter",1]]

具体问题是如何将常规 lat/long 点转换为 shapefile 中的点。

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

感谢任何帮助。

Answer 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()

Answer 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()

将 shapefile 纬度转换为真实纬度

Converting shapefile latitude to real one

python

shapefile