1

I am working with an open data shapefile that uses easting/northing coordinates. My geocoded data set uses lat/long, so I am trying to reproject/transform the shapefile to use lat/long.

```{r}
nwk_shp <- st_read('http://data.ci.newark.nj.us/dataset/61038922-a789-4abe-9928-79f9fe96c307/resource/90ead04f-dffd-4b3a-93c7-39cd4b190d4f/download/superneighborhoodswgs84.geojson') %>% 
  as('Spatial')
```

Reading layer `OGRGeoJSON' from data source `http://data.ci.newark.nj.us/dataset/61038922-a789-4abe-9928-79f9fe96c307/resource/90ead04f-dffd-4b3a-93c7-39cd4b190d4f/download/superneighborhoodswgs84.geojson' using driver `GeoJSON'
Simple feature collection with 21 features and 8 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 561009.6 ymin: 671183.6 xmax: 598979.6 ymax: 712226.2
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs

The shapefile uses northing/easting coordinates:

```{r}
nwk_shp@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
1   0     18         LOWER CLINTON HILL   686419.4  575672.5
2   1      8 MT PLEASANT/LOWER BROADWAY   701184.5  583709.9
3   2     12            LOWER ROSEVILLE   701167.2  577849.7

Here is the projection string on my shapefile:

```{r}
proj4string(nwk_shp)
```
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

Based on this tutorial, and several others, I have tried several variations on spTransform - none of which seem to transform my data into expected lat/longitude coordinates.

For instance:

```{r}

nwk_shp_transformed <- spTransform(nwk_shp, CRS("+proj=longlat +datum=WGS84"))

nwk_shp_transformed@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
 FID MSLINK                       NAME NORTHING_Y EASTING_X
1   0     18         LOWER CLINTON HILL   686419.4  575672.5
2   1      8 MT PLEASANT/LOWER BROADWAY   701184.5  583709.9
3   2     12            LOWER ROSEVILLE   701167.2  577849.7

Depending on what common lat/long CRS I use, I get a 'non-finite transformation detected':

```{r}

nwk_shp_again <- spTransform(nwk_shp, CRS("+init=epsg:4238"))

nwk_shp_again@data[1:3, ] %>% select(-AREA, -PERIMETER, -CATEGORY)
```
35 projected point(s) not finitenon finite transformation detected:
[1] 574389.9 685357.4      Inf      Inf
Error in .spTransform_Polygon(input[[i]], to_args = to_args, from_args = from_args,  : 
  failure in Polygons 1 Polygon 1 points 1

All of this leads me to suspect that there's something about either

  • the proj string of the original shapefile?
  • something omitted from the CRS I am using in my spTransform command?
    at work here.

I've read the 'Overview of Coordinate Reference Systems in R' but am still stuck.

Desired output would be a shapefile / SpatialPolygonsDataFrame with coordinates in lat / longitude - eg, Latitudes around 40 degrees, and Longitudes around -74 degrees.

  • This nwk_shp@data[1:3, ] can't possibly work after reading with sf::st_read. You have a mix of sf and sp code here, can you be sure and choose one or the other - there's a list of analogous functions between the packages here github.com/r-spatial/sf/wiki/Migrating – mdsumner Jan 24 '18 at 22:18
2

It looks like you're correct with your first suspicion, the original proj4string is incorrectly defined as "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0". You can tell because the coordinates don't look like longitude/latitude decimal degrees.

If you knew the correct coordinate system of the original data, you could correct the proj4string like this

new_shp@proj4string <- CRS("+proj=the correct proj parameters")

and then call spTransform like you did

nwk_shp_transformed <- spTransform(nwk_shp, CRS("+proj=longlat +datum=WGS84"))

I tried this, guessing that the correct coordinate system of the original was web mercator, but that didn't seem to be correct either.

Looking at the original data source - http://data.ci.newark.nj.us/dataset/newark-neighborhoods/resource/90ead04f-dffd-4b3a-93c7-39cd4b190d4f - the problem is pretty clearly on their end.

But this might be a working version of the same? http://data.ci.newark.nj.us/dataset/newark-neighborhoods/resource/fcbbb04d-4684-439c-886f-929cb4e7c200

  • This worked wonderfully! For the record, I used proj4string(nwk_shp) <- CRS("+init:epsg:3424") per @mkennedy below – Andrew Jan 25 '18 at 21:53
1

The coordinate values in the nwk_shp appear to be using a New Jersey State Plane zone, units = US survey feet, based on NAD 1983 or a later re-adjustment (like NAD 1983 (2011), which one is impossible to determine by inspection). The EPSG well-known ID for the NAD 1983 version is 3424.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.