Search code examples
pythonlatitude-longitudegeo

Converting geographic coordinates from GEOSTAT to lat and lng

I've found an interesting datasource of European Population which I think could help me in achieving such a map:

enter image description here

The source document GEOSTAT_grid_POP_1K_2011_V2_0_1.csv looks like this:

| 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 |

Geographic coordinates look to be coded in the GRD_ID column this document indicates 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)

I thought I could get lat and long by parsing the strings. For example in 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

but it makes no sense, it does not correspond to a location in Europe !

It may be that it refers to a geographic coordinate system I'm not aware of.

Solution

  • Thanks to Viktor's comment, I found out that the coordinate system used in my file was EPSG:3035

    Based on python's implementation of Proj4, I could achieve a convincing result with the following code:

    #! /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