5

Does anybody know why there a gap appearing between two point layers create with ESRI and QGIS when I a reproject, from the same source file, to Robinson?

The source coordinate system and destination coordinate system share the same datum. It seems there's a shift depending on the technology used to reproject.

Here are the steps to reproduce:

I'm using a shapefile point layer in Lat/Long WGS84 Coordinate System and I reproject it to Robinson (ESRI:54030) using the Tool "Project" in ArcMap

I do the same thing in QGIS. I take the same point layer in WGS84 Lat/Long and reproject it to Robinson (ESRI:54030) using the "Save as.."

Result in QGIS

Download the sample point data here

At the request of mkennedy, here is the cs2cs conversion of regular 10° points:

deg_E   deg_N   sph_X     sph_Y      ell_X       ell_Y
 10       0   943711.34    0        944768.52       0
 10      10   939370.28 1068322.39  940422.59   1069519.16
 10      20   926913.31 2136644.78  927951.66   2139038.32
 10      30   905962.87 3204967.3   906977.76   3208557.61
 10      40   869724.36 4271566.44  870698.65   4276351.58
 10      50   819047.09 5320935.04  819964.61   5326895.73
 10      60   753647.9  6328948.79  754492.16   6336038.68
 10      70   678150.95 7267177.19  678910.63   7275318.12
 10      80   586327.84 8093403.95  586984.66   8102470.44
 10      90   502243.15 8615503.11  502805.77   8625154.47

Scaling the results leads to the factors in https://en.wikipedia.org/wiki/Robinson_projection.

nmtoken
  • 13,355
  • 5
  • 38
  • 87
Jerome Alary
  • 227
  • 1
  • 2
  • 7
  • when I try it directly in GDAL, I get the same results as QGIS. I wonder if ESRI is using a different source projection than 4326, as the degree unit in the prj file is slightly different from that in the 4326 spec (see http://spatialreference.org/ref/epsg/4326/html/)? (0.017453292519943295 vs. 0.01745329251994328) – neuhausr Sep 02 '15 at 14:25
  • How big/long are the differences? – bugmenot123 Sep 02 '15 at 14:46
  • 1
    Distance between points from both layer are roughly from 6 to 8 km – Jerome Alary Sep 02 '15 at 15:04
  • I did try to reproject using the same angular units ((0.017453292519943295 vs. 0.01745329251994328) but i.m getting the same result – Jerome Alary Sep 02 '15 at 15:06
  • The difference in deg-rad conversion isn't enough for mm differences let alone km differences. Let me think about it (I'm the Esri projection specialist). – mkennedy Sep 02 '15 at 16:03
  • 3
    The projections are set up differently--PROJ.4 has pre-calculated values every 5 degrees before it interpolates. Esri uses the table from the original paper, which isn't as dense. I'll have to see what a developer thinks. Esri also has a Robinson (ARC/INFO) version. That gives different answers. It was implemented by a contractor probably 25-30 years ago (before me) and we're still not sure why it's different. – mkennedy Sep 02 '15 at 16:36
  • Can you add the spherical Robinson projection (ESRI 53030) for both software too? – AndreJ Sep 02 '15 at 18:22
  • Yes, I can add EPSG 53030 in both. The Gap result is the same. Same thing with Van der Grinten (EPSG 54029) – Jerome Alary Sep 02 '15 at 18:45
  • @AndreJ but I don't have an answer yet, just speculation! – mkennedy Sep 02 '15 at 20:21
  • Dev thinks the different interpolation grids could be the culprit. I don't have PROJ.4 (or any SW that uses it) running. Do you have any test points at even 10 degree intervals? Those should give the same answers. – mkennedy Sep 02 '15 at 21:02
  • @mkennedy from your answer here http://gis.stackexchange.com/questions/158197/why-do-i-get-correct-area-and-intersect-area-when-use-wrong-projection/158213#158213 I assume that some ESRI ellipsoidal projections are rather spherical ones, with R=a. Can you check if this is true for Robinson too? – AndreJ Sep 03 '15 at 03:50
  • @AndreJ It's true. We thought of that, but afaik, PROJ.4 requires +R*=value, so that shouldn't be a problem. Jerome - the prj/qpj files didn't include the PROJ.4 string. What is it in QGIS? – mkennedy Sep 03 '15 at 13:58
  • Proj.4 definition is +proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs. I'm not sure if they take R=a or R=(a+b)/2 (which is closer to the spherical solution). – AndreJ Sep 03 '15 at 14:50
  • I'm getting much closer answers if I specify a sphere in PROJ.4. I've been trying various values but none are exactly matching the same sphere in ArcGIS. – mkennedy Sep 03 '15 at 17:45

1 Answers1

3

The original paper describing the Robinson projection included a table of values at every 5 degrees along the equator and the central longitude. Those values are used to interpolate intermediate coordinates

This same table is what Esri is using.

The PROJ.4 code has a much expanded table with each 5 degree value having an additional 3 values (possibly done for performance reasons).

An Esri developer thinks that there are different interpolations occurring.

Another possibility is that PROJ.4 is calculating a sphere from the given WGS84 datum while Esri is using the semimajor axis of WGS84, but I haven't been able to check this yet.

Note: I'm getting much closer values (X<1m, Y<120m) if I specify a sphere in PROJ.4 but haven't managed an exact match yet.

mkennedy
  • 18,916
  • 2
  • 38
  • 60
  • Can you add which kind of interpolation ESRI uses between the 5° tabulated values? – AndreJ Sep 07 '15 at 06:13
  • @mkennedy : do you know any way I can match the ESRI Robinson 54030 by using ogr proj4? I need to reproject shapefile using ogr / qgis – Jerome Alary Sep 21 '15 at 19:34