R: Krovak projection in R - identical in ArcGIS, different in R?

Question

I have two files: shapefile and raster, both in the same coordinate system in ArcGIS: EPSG 5514.

Check information on projection in ArcGIS for raster:

S-JTSK_Krovak_East_North
WKID: 5514 Authority: EPSG

Projection: Krovak
False_Easting: 0.0
False_Northing: 0.0
Pseudo_Standard_Parallel_1: 78.5
Scale_Factor: 0.9999
Azimuth: 30.28813975277778
Longitude_Of_Center: 24.83333333333333
Latitude_Of_Center: 49.5
X_Scale: -1.0
Y_Scale: 1.0
XY_Plane_Rotation: 90.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_S_JTSK
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_S_JTSK
  Spheroid: Bessel_1841
    Semimajor Axis: 6377397.155
    Semiminor Axis: 6356078.962818189
    Inverse Flattening: 299.1528128

Check information on projection in ArcGIS for Shapefile:

S-JTSK_Krovak_East_North
WKID: 5514 Authority: EPSG

Projection: Krovak
False_Easting: 0.0
False_Northing: 0.0
Pseudo_Standard_Parallel_1: 78.5
Scale_Factor: 0.9999
Azimuth: 30.28813975277778
Longitude_Of_Center: 24.83333333333333
Latitude_Of_Center: 49.5
X_Scale: -1.0
Y_Scale: 1.0
XY_Plane_Rotation: 90.0
Linear Unit: Meter (1.0)

Geographic Coordinate System: GCS_S_JTSK
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_S_JTSK
  Spheroid: Bessel_1841
    Semimajor Axis: 6377397.155
    Semiminor Axis: 6356078.962818189
    Inverse Flattening: 299.1528128

in ArcGIS totally identical coordinate systems.

---

If I import them in R by

library(sp)
library(rgdal)
library(rgeos)
library(raster)

# import raster
r1<-raster("F:/.../r1_2013.img")

#import shp
d.shp <- readOGR(dsn = "F:/.../dummy_data_arcgis_r",
                     layer = "np_by_manag")

and I want to check data projection, there is difference in projection in my PREVIOUSLY IDENTICAL CRS:

> proj4string(r1)    # RASTER
[1] "+proj=krovak +lat_0=24.83333333333333 +lon_0=49.5 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

> proj4string(d.shp)   # SHAPEFILE
[1] "+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333 +k=0.9999 +x_0=0 +y_0=0 +ellps=bessel +units=m +no_defs"
>

which cause problems in any spatial operation I want to perform further.

Why R modify my projection? Is it caused by various nature of my files: shapefile / raster? How can I get rid of it and which one of Krovak CRS is correct?

My aim, after finishing this transformations, is to:

1. convert raster r1 to polygon p.r1,
2. intersect polygon p.r1 with d.shp
3. calculate areas of each polygon

---

# convert raster to polygon
p.r1<-rasterToPolygons(r1, dissolve = T)

# intersect polygon from raster with another one
int.r <-intersect(p.r1,d.shp)  

# calculate areas (in m2) for every polygon in intersected shp
int.r$area <- sapply(int.r@polygons, function(x) x@Polygons[[1]]@area)

My dummy files are attached here: http://ulozto.sk/xKoewAzg/dummy-data-arcgis-r-zip

Answer 1

Looks like ESRI decided to define their own projection parameters:

http://osgeo-org.1560.x6.nabble.com/Question-about-Krovak-projection-and-ESRI-XY-Plane-Rotation-parameter-td4291471.html

There has also been some discussion on the GDAL list recently:

http://comments.gmane.org/gmane.comp.gis.gdal.devel/42132

I suspect the raster package is applying some of the extra parameters in the projection, and possibly doing it wrong.

The underlying extents in the data are the same:

> extent(r1)
class       : Extent 
xmin        : -373118.6 
xmax        : -366218.6 
ymin        : -1176124 
ymax        : -1170154 
> extent(d.shp)
class       : Extent 
xmin        : -373118.6 
xmax        : -366212.9 
ymin        : -1176124 
ymax        : -1170145 

and if you transform the shape back to lat-long you get something in the Czech and Slovakian region:

> bbox(spTransform(d.shp,CRS("+init=epsg:4326")))
       min      max
x 19.70567 19.79555
y 49.23434 49.29116

EXCEPT that the x and y coordinates are swapped! Prague is about 14N 50E. This appears to be the major problem with this CRS - the coordinates are in the wrong order and one is increasing in a southerly direction.

The projection read in by raster is even wronger, when projected to lat-long the coordinates are in completely the wrong place:

> extent(projectExtent(r1,"+init=epsg:4326"))
class       : Extent 
xmin        : 44.3646 
xmax        : 44.46467 
ymin        : 49.24811 
ymax        : 49.30581 

I think if you set the raster projection to the same as the shape, then you can do your spatial overlays etc because nothing will complain about mismatching coordinate systems:

> projection(r1) = projection(d.shp)
> 

But bear in mind that what you plot on the screen is possibly rotated and reflected compared to the ground. However if I load them into QGIS and give them both an EPSG of 5514, then when compared to how R plots them, I get this:

R on the left, QGIS on the right, where your data appears to be in a forest on the Poland/Slovakia border. QGIS has a definition for EPSG:5514, unlike R. The QGIS map looks the same with EPSG:5514, EPSG:2065 or EPSG:5221, the last two of which I do have in R, so I could assign those to my data. I don't see much difference between them, possibly a tiny survey shift that doesn't materially affect the data.

So the raster is definitely wrong. In fact when I do what I should have done at the start and looked at the projection parameters of the file, it DOESNT have the same projection parameters as the shape:

$ gdalinfo r1_2013.img 
Driver: HFA/Erdas Imagine Images (.img)
Files: r1_2013.img
       r1_2013.img.aux.xml
Size is 230, 199
Coordinate System is:
PROJCS["S-JTSK_Krovak_East_North",
    GEOGCS["S_JTSK",
        DATUM["S_JTSK",
            SPHEROID["Bessel_1841",6377397.155,299.1528128000033],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Krovak"],
    PARAMETER["latitude_of_center",24.83333333333333],
    PARAMETER["longitude_of_center",49.5],
    PARAMETER["azimuth",30.28813975277778],

Notice it has the central longitude at 49 degrees rather than 24.83. How did that happen? Did ArcGIS write the wrong CRS to the image file??

In short, as long as the back-transform to lat-long looks sensible, it probably doesn't matter.