将地理坐标从 GEOSTAT 转换为纬度和经度

Question

我发现了一个有趣的 datasource of European Population，我认为它可以帮助我实现这样的地图：

源文档 GEOSTAT_grid_POP_1K_2011_V2_0_1.csv 如下所示：

| TOT_P | GRD_ID | CNTR_CODE | METHD_CL | YEAR | DATA_SRC | TOT_P_CON_DT | |-------|---------------|-----------|----------|------|----------|--------------| | 8 | 1kmN2689E4337 | DE | A | 2011 | DE | other | | 7 | 1kmN2689E4341 | DE | A | 2011 | DE | other |

地理坐标看起来编码在本文档指示的 GRD_ID 列中 Appendix1_WP1C_production-procedures-bottom-up.pdf:

Grid cell identification codes are based on grid cell’s lower left-hand corner coordinates truncated by grid cell size (e.g. 1kmN4534E5066 is result from coordinates Y=4534672, X=5066332 and the cell size 1000)

我以为我可以通过解析字符串来获取纬度和经度。例如 Python:

import re
string = "1kmN2691E4341"
lat = float(re.sub('.*N([0-9]+)[EW].*', '\1', string))/100
lng = float(re.sub('.*[EW]([0-9]+)', '\1', string))/100

print lat, ",", lng

Output 26.91 , 43.41

但是没有意义，它不对应欧洲的位置！

可能是指我不知道的地理坐标系。

Answer 1

感谢 Viktor 的评论，我发现我的文件中使用的坐标系是 EPSG:3035

基于python对Proj4的实现，我可以用下面的代码得到一个令人信服的结果：

#! /usr/bin/python
# coding: utf-8

import re
from pyproj import Proj, transform

string = "1kmN2326E3989"

x1 = int(re.sub('.*[EW]([0-9]+)', '\1', string))*1000
y1 = int(re.sub('.*N([0-9]+)[EW].*', '\1', string))*1000

inProj = Proj(init='EPSG:3035')
outProj = Proj(init='epsg:4326')

lng,lat = transform(inProj,outProj,x1,y1)
print lat,lng

Output : 43.9613760836 5.870517281

将地理坐标从 GEOSTAT 转换为纬度和经度

Converting geographic coordinates from GEOSTAT to lat and lng

python

geo

latitude-longitude