0

I use data from NOAA for some analysis in R and I want to transform coordinates to EPSG:54004 or something useful.

There is something what they write about their coordinates.

...    
POS: 29-34 
 GEOPHYSICAL-POINT-OBSERVATION latitude coordinate 
 The latitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where southern 
 hemisphere is negative. 
 MIN: -90000 MAX: +90000 
 UNITS: Angular Degrees 
 SCALING FACTOR: 1000 
 DOM: A general domain comprised of the numeric characters (0-9), a plus 
 sign (+), and a minus sign (-). 
 +99999 = Missing 

POS: 35-41 
 GEOPHYSICAL-POINT-OBSERVATION longitude coordinate 
 The longitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where values west from 
 000000 to 179999 are signed negative. 
 MIN: -179999 MAX: +180000 UNITS: Angular Degrees 
 SCALING FACTOR: 1000 
 DOM: A general domain comprised of the numeric characters (0-9), a plus 
 sign (+), and a minus sign (-). 
 +999999 = Missing
... 

And my problem is, that I can't transform this coodrinates correctly. I use this R command:

    sc <- cbind(st$LAT, st$LON)
    ptransform(sc/180*pi, 
            '+proj=latlong +ellps=sphere',
            '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

But coordinates are somehow wrong as you can see. The points should be in borders because they are meteostations from this countries. Map with points.

  • Which coordinate system (EPSG code) would you like to transform to? If you are familiar with python, I can post a script that does that in a few lines of code – Matej May 07 '14 at 00:22

2 Answers2

1

I solved the problem, it was in input to ptransform function in R.
Instead sc <- cbind(st$LAT, st$LON) I have to use sc <- cbind(st$LON, st$LAT). Then as you say I have to use EPSG:4326 as input CRS and EPSG:3857 as destination CRS. And the transformation with this command.

    tr <- ptransform(sc/180*pi, 
            '+proj=longlat',
            '+proj=merc') 

Correct transformation

0

In Python using pyproj, this is very simple:

import pyproj

def transform(epsg_in, epsg_out, x_in, y_in):
    # define source and destination coordinate systems based on the ESPG code
    srcProj = pyproj.Proj(init='epsg:%i' % epsg_in, preserve_units=True)
    dstProj = pyproj.Proj(init='epsg:%i' % epsg_out, preserve_units=True)

    # perform the transformation
    x_out,y_out = pyproj.transform(srcProj, dstProj, x_in, y_in)
    print '(%.5f, %.5f) EPSG: %i => (%.5f, %.5f) EPSG: %i' % (x_in, y_in, epsg_in, x_out, y_out, epsg_out)
    return x_out, y_out

transform(**{
    "epsg_in" : 4326,
    "epsg_out" : 3912, 
    "x_in" : 14.5,
    "y_in" : 46.05
})

You only need to know epsg codes and input coordinates.

Matej
  • 1,758
  • 1
  • 13
  • 20