0

I'm trying to use st_transform (from the R package sf) on a sf object and am getting the following error when projecting from WGS 1984 to Albers Equal Area:

> sf_obj_prj <- sf_obj %>% 
   st_transform(crs = 
     "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")
OGR: Not enough data 
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
GDAL Error 1: Length of input WKB is too small

I have checked that my map has a CRS:

> st_crs(sf_obj)
Coordinate Reference System:
EPSG: 4326 
proj4string: "+proj=longlat +datum=WGS84 +no_defs"

I have also checked that all polygons are valid:

> which(!is_true(st_is_valid(sf_obj)))
integer(0)

(is_true is my user-defined function that returns TRUE if true and FALSE if NA/false)

I've tried searching but could not find anything about this error.

EDIT 1: It seems that this is arising after using rbind to bind two objects.

Here is the st_bbox on the two objects:

> st_bbox(sf1)
      xmin       ymin       xmax       ymax 
-73.991244 -11.145219 -66.619911  -7.111824 
> st_bbox(sf2)
     xmin      ymin      xmax      ymax 
-62.68262  80.03352 -57.99074  82.14633 

st_transform works on the objects separately, but not when I use rbind:

> x = st_transform(sf1, crs = prj_alber)
> y = st_transform(sf2, crs = prj_alber)
> z = st_transform(rbind(sf1, sf2), crs = prj_alber)
OGR: Not enough data 
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
  GDAL Error 1: Length of input WKB is too small

EDIT 2: Here is code to reproduce the error (note: requires downloading data from the Brazilian statistics office)

download.file(
  "http://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2000/setor_rural/projecao_geografica/censo_2000/e1000_generico_dxf/uf/ac.zip",
  paste0(directory, "sf1.zip"))
unzip(paste0(directory, "sf1.zip"), exdir = paste0(directory, "sf1"))

download.file(
  "http://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2000/setor_urbano/ac/1200203/1200203.zip",
  paste0(directory, "sf2.zip"))
unzip(paste0(directory, "sf2.zip"), exdir = paste0(directory, "sf2"))


sf1 <- st_read(paste0(directory, "sf1/AC/12SE1000G.dxf"))
# (original crs undefined)
st_crs(sf1) <- "+init=epsg:5527"
sf1 <- st_transform(sf1, crs = "+init=epsg:4326")
sf1 <- sf1 %>% dplyr::select(ID = Text)

sf2 <- st_read(paste0(directory, "sf2/1200203.SHP"))
# (original crs was defined for northern hemisphere, this is the epsg code for southern hemisphere)
st_crs(sf2) <- "+init=epsg:32718"
sf2 <- sf2 %>% st_transform(sf2, crs = "+init=epsg:4326")
# (some polygons were invalid so fixing them here)
invalid = which(!is_true(st_is_valid(sf2)))
sf2[invalid, ] = st_make_valid(sf2[invalid,])
sf2 <- sf2 %>% dplyr::select(ID = ID_)

> x = st_transform(sf1, crs = prj_alber)
> y = st_transform(sf2, crs = prj_alber)
> z = st_transform(rbind(sf1, sf2), crs = prj_alber)
OGR: Not enough data 
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
  GDAL Error 1: Length of input WKB is too small

But as per Spacedman's solution below, removing the extra dimension fixes the issue:

z <- st_transform(st_zm(rbind(sf1, sf2)), crs = prj_alber)
  • Works for me for a test epsg 4326 dataset nc derived from running example(st_read). Must be something in your data. Are the source coordinates in the right range? What's st_bbox telling you? Does it work on a subset of features in your source data? – Spacedman Mar 24 at 10:27
  • Thanks for the questions - I've updated the question to elaborate more about the data I'm working with – mochi Mar 24 at 19:16
  • Very very weird. Any chance you can share your data? Its not a general problem since I just rbinded two halves of an sf data frame and it works fine. Does it fail with other projections (try "+init=epsg:3857")? Can you show summary() of your objects if you cant share them – Spacedman Mar 24 at 19:24
  • I can share - it's publicly available data from the Brazilian government. I'm not sure how to do it though, how would you suggest I share? – mochi Mar 24 at 19:37
  • Sorry - I am being silly. Here are the links to the datasets. This is sf1: geoftp.ibge.gov.br/organizacao_do_territorio/… and this is sf2: geoftp.ibge.gov.br/organizacao_do_territorio/… – mochi Mar 24 at 19:40
0

I can get this error if I mix objects with different dimensions. Here's a reproducible example:

First make a data frame with two points:

d = st_as_sf(data.frame(x=1:2,y=1:2,z=c("A","B")), coords=1:2)

Next make two points with different dimensionalities:

p1 = st_point(c(1,1,2,3), dim="XYZM")
p2 = st_point(c(1,2), dim="XY")

Replace the points in the data frame (at which point I think there's probably an easier way to do this by nm):

d$geometry[1] = p1
d$geometry[2] = p2

Give the thing a CRS - lat-long is valid:

st_crs(d)=4326

And now transform to oooh, Google Web Mercator:

st_transform(d,3857)
# OGR: Not enough data 
# Error in CPL_transform(x, crs$proj4string) : OGR error

To fix, wrap the object in st_zm which removes any Zs or Ms and returns a purely XY object:

> st_transform(st_zm(d),3857)
Simple feature collection with 2 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 111319.5 ymin: 111325.1 xmax: 111319.5 ymax: 222684.2
epsg (SRID):    3857
proj4string:    +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs
  z                  geometry
1 A POINT (111319.5 111325.1)
2 B POINT (111319.5 222684.2)

If you need the Zs and the Ms then you'll have to take them from the original objects and add them back. Remembering that some of the rows didn't have them in the first place....

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.