Re: [R-sig-Geo] Adding spatial tables to existing SpatiaLite DBs

[Roger Bivand]

On Tue, 3 Nov 2015, Barry Rowlingson wrote:


I don't think what I'm trying to do is "appending". I'm trying to
write two spatial data tables with different names in the single
spatial database file. The database file is the DSN and the tables are
layers - and rgdal is quite happy to add Shapefile "layers" to a
folder "DSN"...

Have just upgraded to rgdal 1.1-1 and still the same problem. Here's a
reprex without prompts:

require(rgdal)
pts=data.frame(x=runif(10),y=runif(10),z=1:10)
coordinates(pts)=~x+y
file.remove("tmpfile.db")

# write layer pts, this works
writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

# write layer pts2
writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

Interestingly when I try and write with the *same* table name, rgdal
helpfully suggests "layer exists, use a new layer name", but when I
obey, I get the error behaviour I've described


There is nothing (much) in rgdal::writeOGR() that knows much about the 
drivers and (see separate recent thread) the implementation is really 
based on the simplest formats that existed 10 years ago. Spatialite is 
fast-changing, and not very stable in my experience. You're asking for OGR 
to handle this, but writeOGR() isn't aware that it needs to ask OGR for 
things that were not there when it was written (it was based on GRASS 
v.out.ogr from even earlier). It would be nice if it worked by chance, 
unspecified, but more needs to be done to check, and it looks as though 
it doesn't work (the single added layer overwrites the existing layer, 
problably because the dsn= isn't opened to append).


Does the same thing happen without specifying Spatialite? Does the same 
thing happen with PostGIS? Other DBs? We know that when dsn= is a 
directory and the driver is ESRI Shapefile, it does what we expect, but 
should we expect it to do that when dsn= is a file?


Please look around line 60 in rgdal/R/ogr_write.R to see how writeOGR() 
handles dsn and layer checking. This may need conditioning on the driver - 
there is already a kludgy "fix" for dsn= deletion for shapefiles for GDAL 

=2.


When you're checking improvements to writeOGR(), please set a baseline 
using gdalUtils::ogr2ogr() so that we know where we are. We'll need to 
support GDAL < 2 and GDAL >= 2, which use drivers differently.


Best wishes,

Roger



Barry



On Mon, Nov 2, 2015 at 11:21 PM, Roger Bivand  wrote:

On Tue, 3 Nov 2015, Barry Rowlingson wrote:


I can create a SpatiaLite DB file and put a layer in it, but if I try
and add another layer, rgdal fails. Example:

Versions etc:


require(rgdal)

Loading required package: rgdal
Loading required package: sp
prgdal: version: 1.0-7, (SVN revision 559)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
Path to GDAL shared files: /usr/share/gdal/1.11
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)
Linking to sp version: 1.2-0

Create a simple points data set, write it:


pts = data.frame(x=runif(10),y=runif(10),z=1:10)
coordinates(pts)=~x+y
writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",

dataset_option="SPATIALITE=YES")

Note the use of the dataset_option to make this a proper SpatiaLite,
and not just an SQLite table. The output file is about 4Mb and has a
lot of metadata tables in it. The file loads into QGIS which
recognises it as a SpatiaLite table and I can plot the points.

Now try and create another spatial table (this time, "pts2")  in the
same database file:


writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",

dataset_option="SPATIALITE=YES")
Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
dataset_option = "SPATIALITE=YES") :
 Creation of output file failed


No reported error in 1.1-1, but the original layer is overwritten. The
layer in the existing dsn is overwritten. Appending is not something
writeOGR knows about:


writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",

dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",

dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
dataset_option = "SPATIALITE=YES",  :
   layer exists, use a new layer name

ogrListLayers("tmpfile.db")

[1] "pts2"
attr(,"driver")
[1] "SQLite"
attr(,"nlayers")
[1] 1

For now, system("ogr2ogr ..."), or something from gdalUtils? Contribution
to writeOGR() to implement appending as in
http://www.gdal.org/drv_sqlite.html?

Roger



I've tried various other options with no success. Is this possible?
Or do I ditch writeOGR for this and create the table by hand... ick...

Barry

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org

Re: [R-sig-Geo] Adding spatial tables to existing SpatiaLite DBs

[Roger Bivand]

On Tue, 3 Nov 2015, Roger Bivand wrote:


On Tue, 3 Nov 2015, Barry Rowlingson wrote:


 I don't think what I'm trying to do is "appending". I'm trying to
 write two spatial data tables with different names in the single
 spatial database file. The database file is the DSN and the tables are
 layers - and rgdal is quite happy to add Shapefile "layers" to a
 folder "DSN"...

 Have just upgraded to rgdal 1.1-1 and still the same problem. Here's a
 reprex without prompts:

 require(rgdal)
 pts=data.frame(x=runif(10),y=runif(10),z=1:10)
 coordinates(pts)=~x+y
 file.remove("tmpfile.db")

 # write layer pts, this works
 writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
 dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

 # write layer pts2
 writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
 dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

 Interestingly when I try and write with the *same* table name, rgdal
 helpfully suggests "layer exists, use a new layer name", but when I
 obey, I get the error behaviour I've described


There is nothing (much) in rgdal::writeOGR() that knows much about the 
drivers and (see separate recent thread) the implementation is really based 
on the simplest formats that existed 10 years ago. Spatialite is 
fast-changing, and not very stable in my experience. You're asking for OGR to 
handle this, but writeOGR() isn't aware that it needs to ask OGR for things 
that were not there when it was written (it was based on GRASS v.out.ogr from 
even earlier). It would be nice if it worked by chance, unspecified, but more 
needs to be done to check, and it looks as though it doesn't work (the single 
added layer overwrites the existing layer, problably because the dsn= isn't 
opened to append).


Does the same thing happen without specifying Spatialite? Does the same thing 
happen with PostGIS? Other DBs? We know that when dsn= is a directory and the 
driver is ESRI Shapefile, it does what we expect, but should we expect it to 
do that when dsn= is a file?


Please look around line 60 in rgdal/R/ogr_write.R to see how writeOGR() 
handles dsn and layer checking. This may need conditioning on the driver - 
there is already a kludgy "fix" for dsn= deletion for shapefiles for GDAL 

=2.


When you're checking improvements to writeOGR(), please set a baseline using 
gdalUtils::ogr2ogr() so that we know where we are. We'll need to support GDAL 
< 2 and GDAL >= 2, which use drivers differently.


tf1 <- tempfile()
tf2 <- tempfile()
require(rgdal)
pts=data.frame(x=runif(10),y=runif(10),z=1:10)
coordinates(pts)=~x+y
writeOGR(pts, tf1, "pts", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
writeOGR(pts, tf2, "pts2", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
ogrListLayers(tf1)
ogrListLayers(tf2)
library(gdalUtils)
ogr2ogr(tf2, tf1, append=TRUE, nln="pts2a")
ogrListLayers(tf1)
ogrInfo(tf1, "pts")
ogrInfo(tf1, "pts2a")

Roger



Best wishes,

Roger



 Barry



 On Mon, Nov 2, 2015 at 11:21 PM, Roger Bivand  wrote:
>  On Tue, 3 Nov 2015, Barry Rowlingson wrote:
> 
> >  I can create a SpatiaLite DB file and put a layer in it, but if I try

> >  and add another layer, rgdal fails. Example:
> > 
> >  Versions etc:
> > 
> > >  require(rgdal)

> >  Loading required package: rgdal
> >  Loading required package: sp
> >  prgdal: version: 1.0-7, (SVN revision 559)
> >  Geospatial Data Abstraction Library extensions to R successfully 
> >  loaded

> >  Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
> >  Path to GDAL shared files: /usr/share/gdal/1.11
> >  Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
> >  Path to PROJ.4 shared files: (autodetected)
> >  Linking to sp version: 1.2-0
> > 
> >  Create a simple points data set, write it:
> > 
> > >  pts = data.frame(x=runif(10),y=runif(10),z=1:10)

> > >  coordinates(pts)=~x+y
> > >  writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
> >  dataset_option="SPATIALITE=YES")
> > 
> >  Note the use of the dataset_option to make this a proper SpatiaLite,

> >  and not just an SQLite table. The output file is about 4Mb and has a
> >  lot of metadata tables in it. The file loads into QGIS which
> >  recognises it as a SpatiaLite table and I can plot the points.
> > 
> >  Now try and create another spatial table (this time, "pts2")  in the

> >  same database file:
> > 
> > >  writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",

> >  dataset_option="SPATIALITE=YES")
> >  Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
> >  dataset_option = "SPATIALITE=YES") :
> >   Creation of output file failed
> 
>  No reported error in 1.1-1, but the original layer is overwritten. The

>  layer in the existing dsn is overwritten. Appending is not something
>  writeOGR knows about:
> 
> >  writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",

>  dataset_option="SPATIALITE=YES", 

Re: [R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR

[Roger Bivand]

Hi Andy,

On Mon, 26 Oct 2015, Andy Teucher wrote:


Hi Roger,

Sorry I haven't been much help here - it's not for lack of trying,
just lack of skills. I've been poring through OGR_write.cpp and trying
to figure out language and the flow of things... a few thoughts
occurred to me (which are probably obvious to you):



Your example was very helpful, the problem was that when writeOGR() and 
the underlying compiled code was written, the model used was that of the 
ESRI Shapefile driver, rather than OGC SFS, which was first published in 
2004. Most older formats (and their drivers in OGR) use a OGR geometry 
factory function to organize polygons on reading - but this is not present 
in shapelib (bundled with the R maptools package), so when designing 
polygon classes for sp, and writing rgdal::writeOGR() about 10 years ago, 
this was not known, and was not revisited subsequently, until you happily 
spotted the bug.


GEOS in rgeos uses OGC SFS definitions (the ones you are using below), but 
the SpatialPolygons/Polygons/Polygon structure was already established by 
then, and a fix (using comment attributes in Polygons objects) was used, 
rather than a major change to Polygons objects (in 2010).


The revision to rgdal::writeOGR() now branches on the input objects - if 
they have comments (applied in code, not user-facing), these are trusted, 
and the MultiPolygon objects constructed containing one or more Polygon 
objects with one exterior ring and zero or more interior rings on that 
basis. However, if the input objects do not have comment attributes, 
OGRGeometryFactory::organizePolygons() (C++) is used to create properly 
organized MultiPolygon objects for export.


This means that the pre-OGC SFS drivers continue to work, but that the OGC 
SFS drivers, like GeoJSON, now also work correctly when the second ring in 
an object is not assumed to be an interior ring (rgdal < 1.1-1). A test 
set for all the possible cases has been introduced.


The revision is now in source package form on CRAN:

https://cran.r-project.org/web/packages/rgdal/index.html

and a Windows binary package will be available later on Tuesday I hope. 
The OSX revision may take a little longer, but we are grateful to have 
what we have anyway.


Thanks for a clear and very helpful bug report.

Best wishes,

Roger


1) Is this bug because holes are standalone polygons in sp objects,
while they are members of polygon objects in geojson, wkt, etc?

2) On a related note, I notice that a single polygon with a hole is
identified as a multipolygon containing a single polygon with a hole,
when using writeOGR and one of the affected drivers.  Practically I
don't think this matters too much, but it's not quite right.

3) In lines 116-160 the polygon/multipolygon and line/multiline
determination is made for the entire layer. Can/should this be made
for each feature in the layer, as you could have a layer with a
mixture of single polygons and multipolygons?

4) I feel like there ought to be one more layer of looping through the
polyons. e.g.:

loop over each feature;
determine whether each feature is a single polygon (length 1 or
all rings but one are holes), or a multipolygon:
   if a polygon loop over and add ring(s);
   if a multipolygon loop over polygons:
   for each polygon add ring(s) and somehow assign holes to
parent polygons

Again, apologies if I'm being totally obvious (or way off base).
Thanks for all your work on these packages, and on this issue.

Andy


On Thu, Oct 22, 2015 at 11:07 PM, Roger Bivand  wrote:

On Fri, 23 Oct 2015, Andy Teucher wrote:


On my Mac, I am running rgdal 1.0-7 and GDAL 1.11.3


sessionInfo()



OK, thanks.

The problem also affects the SQLite driver (the CRAN binaries do not have
this driver). Could someone please check the PostGIS driver - my guess is
that all OGC SFS compliant drivers may be affected.

Could somebody also please check whether this is a recently introduced issue
or not? For example somebody still on rgdal < 1.0? src/OGR_write.cpp was
changed 2015-08-21, but the last previous change was 2015-06-11, then
2015-05-31, 2014-08-17 ... from the ChangeLog visible on the package CRAN
page.

Roger







R version 3.2.2 (2015-08-14)

Platform: x86_64-apple-darwin13.4.0 (64-bit)

Running under: OS X 10.11 (El Capitan)




locale:

[1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8




attached base packages:

[1] stats graphics  grDevices utils datasets  methods   base




other attached packages:

[1] rgdal_1.0-7 sp_1.2-1




loaded via a namespace (and not attached):

[1] tools_3.2.2 grid_3.2.2  lattice_0.20-33




Andy

On Thu, Oct 22, 2015 at 11:28 AM, Andy Teucher 
wrote:


One more test: WKT has the same issue:
library(rgdal)
js <- '{
"type": "MultiPolygon",
"coordinates": 102.0, 2.0], [103.0, 2.0], [103.0, 3.0], [102.0, 3.0],
[102.0, 2.0]]], [[[100.0, 0.0], [101.0, 0.0], [101.0, 1.0], [100.0, 

Re: [R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR

[Andy Teucher]

Amazing (and fast) work Roger, and thanks for the great explanation.
I can see why you wouldn't want to redo the way Polygons objects are
structured - I'm sure it would be break a huge amount of existing
code.

Just curious, does OSX take longer because you have to deal with
multiple possible GDAL versions?

Thanks again,
Andy

On Mon, Nov 2, 2015 at 1:21 PM, Roger Bivand  wrote:
> Hi Andy,
>
> On Mon, 26 Oct 2015, Andy Teucher wrote:
>
>> Hi Roger,
>>
>> Sorry I haven't been much help here - it's not for lack of trying,
>> just lack of skills. I've been poring through OGR_write.cpp and trying
>> to figure out language and the flow of things... a few thoughts
>> occurred to me (which are probably obvious to you):
>>
>
> Your example was very helpful, the problem was that when writeOGR() and the
> underlying compiled code was written, the model used was that of the ESRI
> Shapefile driver, rather than OGC SFS, which was first published in 2004.
> Most older formats (and their drivers in OGR) use a OGR geometry factory
> function to organize polygons on reading - but this is not present in
> shapelib (bundled with the R maptools package), so when designing polygon
> classes for sp, and writing rgdal::writeOGR() about 10 years ago, this was
> not known, and was not revisited subsequently, until you happily spotted the
> bug.
>
> GEOS in rgeos uses OGC SFS definitions (the ones you are using below), but
> the SpatialPolygons/Polygons/Polygon structure was already established by
> then, and a fix (using comment attributes in Polygons objects) was used,
> rather than a major change to Polygons objects (in 2010).
>
> The revision to rgdal::writeOGR() now branches on the input objects - if
> they have comments (applied in code, not user-facing), these are trusted,
> and the MultiPolygon objects constructed containing one or more Polygon
> objects with one exterior ring and zero or more interior rings on that
> basis. However, if the input objects do not have comment attributes,
> OGRGeometryFactory::organizePolygons() (C++) is used to create properly
> organized MultiPolygon objects for export.
>
> This means that the pre-OGC SFS drivers continue to work, but that the OGC
> SFS drivers, like GeoJSON, now also work correctly when the second ring in
> an object is not assumed to be an interior ring (rgdal < 1.1-1). A test set
> for all the possible cases has been introduced.
>
> The revision is now in source package form on CRAN:
>
> https://cran.r-project.org/web/packages/rgdal/index.html
>
> and a Windows binary package will be available later on Tuesday I hope. The
> OSX revision may take a little longer, but we are grateful to have what we
> have anyway.
>
> Thanks for a clear and very helpful bug report.
>
> Best wishes,
>
> Roger
>
>
>> 1) Is this bug because holes are standalone polygons in sp objects,
>> while they are members of polygon objects in geojson, wkt, etc?
>>
>> 2) On a related note, I notice that a single polygon with a hole is
>> identified as a multipolygon containing a single polygon with a hole,
>> when using writeOGR and one of the affected drivers.  Practically I
>> don't think this matters too much, but it's not quite right.
>>
>> 3) In lines 116-160 the polygon/multipolygon and line/multiline
>> determination is made for the entire layer. Can/should this be made
>> for each feature in the layer, as you could have a layer with a
>> mixture of single polygons and multipolygons?
>>
>> 4) I feel like there ought to be one more layer of looping through the
>> polyons. e.g.:
>>
>> loop over each feature;
>> determine whether each feature is a single polygon (length 1 or
>> all rings but one are holes), or a multipolygon:
>>if a polygon loop over and add ring(s);
>>if a multipolygon loop over polygons:
>>for each polygon add ring(s) and somehow assign holes to
>> parent polygons
>>
>> Again, apologies if I'm being totally obvious (or way off base).
>> Thanks for all your work on these packages, and on this issue.
>>
>> Andy
>>
>>
>> On Thu, Oct 22, 2015 at 11:07 PM, Roger Bivand 
>> wrote:
>>>
>>> On Fri, 23 Oct 2015, Andy Teucher wrote:
>>>
 On my Mac, I am running rgdal 1.0-7 and GDAL 1.11.3


 sessionInfo()

>>>
>>> OK, thanks.
>>>
>>> The problem also affects the SQLite driver (the CRAN binaries do not have
>>> this driver). Could someone please check the PostGIS driver - my guess is
>>> that all OGC SFS compliant drivers may be affected.
>>>
>>> Could somebody also please check whether this is a recently introduced
>>> issue
>>> or not? For example somebody still on rgdal < 1.0? src/OGR_write.cpp was
>>> changed 2015-08-21, but the last previous change was 2015-06-11, then
>>> 2015-05-31, 2014-08-17 ... from the ChangeLog visible on the package CRAN
>>> page.
>>>
>>> Roger
>>>
>>>




 R version 3.2.2 (2015-08-14)

 Platform: x86_64-apple-darwin13.4.0 (64-bit)

 

[R-sig-Geo] Adding spatial tables to existing SpatiaLite DBs

[Barry Rowlingson]

I can create a SpatiaLite DB file and put a layer in it, but if I try
and add another layer, rgdal fails. Example:

Versions etc:

 > require(rgdal)
Loading required package: rgdal
Loading required package: sp
prgdal: version: 1.0-7, (SVN revision 559)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
 Path to GDAL shared files: /usr/share/gdal/1.11
 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-0

Create a simple points data set, write it:

 > pts = data.frame(x=runif(10),y=runif(10),z=1:10)
 > coordinates(pts)=~x+y
 > writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
dataset_option="SPATIALITE=YES")

Note the use of the dataset_option to make this a proper SpatiaLite,
and not just an SQLite table. The output file is about 4Mb and has a
lot of metadata tables in it. The file loads into QGIS which
recognises it as a SpatiaLite table and I can plot the points.

Now try and create another spatial table (this time, "pts2")  in the
same database file:

 > writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
dataset_option="SPATIALITE=YES")
Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
dataset_option = "SPATIALITE=YES") :
  Creation of output file failed

 I've tried various other options with no success. Is this possible?
Or do I ditch writeOGR for this and create the table by hand... ick...

Barry

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Adding spatial tables to existing SpatiaLite DBs

[Barry Rowlingson]

I don't think what I'm trying to do is "appending". I'm trying to
write two spatial data tables with different names in the single
spatial database file. The database file is the DSN and the tables are
layers - and rgdal is quite happy to add Shapefile "layers" to a
folder "DSN"...

Have just upgraded to rgdal 1.1-1 and still the same problem. Here's a
reprex without prompts:

require(rgdal)
pts=data.frame(x=runif(10),y=runif(10),z=1:10)
coordinates(pts)=~x+y
file.remove("tmpfile.db")

# write layer pts, this works
writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

# write layer pts2
writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")

Interestingly when I try and write with the *same* table name, rgdal
helpfully suggests "layer exists, use a new layer name", but when I
obey, I get the error behaviour I've described

Barry



On Mon, Nov 2, 2015 at 11:21 PM, Roger Bivand  wrote:
> On Tue, 3 Nov 2015, Barry Rowlingson wrote:
>
>> I can create a SpatiaLite DB file and put a layer in it, but if I try
>> and add another layer, rgdal fails. Example:
>>
>> Versions etc:
>>
>> > require(rgdal)
>> Loading required package: rgdal
>> Loading required package: sp
>> prgdal: version: 1.0-7, (SVN revision 559)
>> Geospatial Data Abstraction Library extensions to R successfully loaded
>> Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
>> Path to GDAL shared files: /usr/share/gdal/1.11
>> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
>> Path to PROJ.4 shared files: (autodetected)
>> Linking to sp version: 1.2-0
>>
>> Create a simple points data set, write it:
>>
>> > pts = data.frame(x=runif(10),y=runif(10),z=1:10)
>> > coordinates(pts)=~x+y
>> > writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
>> dataset_option="SPATIALITE=YES")
>>
>> Note the use of the dataset_option to make this a proper SpatiaLite,
>> and not just an SQLite table. The output file is about 4Mb and has a
>> lot of metadata tables in it. The file loads into QGIS which
>> recognises it as a SpatiaLite table and I can plot the points.
>>
>> Now try and create another spatial table (this time, "pts2")  in the
>> same database file:
>>
>> > writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
>> dataset_option="SPATIALITE=YES")
>> Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
>> dataset_option = "SPATIALITE=YES") :
>>  Creation of output file failed
>
> No reported error in 1.1-1, but the original layer is overwritten. The
> layer in the existing dsn is overwritten. Appending is not something
> writeOGR knows about:
>
>> writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
> dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
>> writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
> dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
> Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
> dataset_option = "SPATIALITE=YES",  :
>layer exists, use a new layer name
>> ogrListLayers("tmpfile.db")
> [1] "pts2"
> attr(,"driver")
> [1] "SQLite"
> attr(,"nlayers")
> [1] 1
>
> For now, system("ogr2ogr ..."), or something from gdalUtils? Contribution
> to writeOGR() to implement appending as in
> http://www.gdal.org/drv_sqlite.html?
>
> Roger
>
>>
>> I've tried various other options with no success. Is this possible?
>> Or do I ditch writeOGR for this and create the table by hand... ick...
>>
>> Barry
>>
>> ___
>> R-sig-Geo mailing list
>> R-sig-Geo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Department of Economics, Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 91 00
> e-mail: roger.biv...@nhh.no
>

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR

[Forrest Stevens]

This is fantastic work, thank you both for troubleshooting it and the fix
Roger. Your tireless work is greatly appreciated! I just ran across this
problem a couple of weeks ago when exporting from R to GeoJSON, and will
happily test out the new version.

Sincerely,
Forrest

On Mon, Nov 2, 2015 at 5:44 PM Andy Teucher  wrote:

> Amazing (and fast) work Roger, and thanks for the great explanation.
> I can see why you wouldn't want to redo the way Polygons objects are
> structured - I'm sure it would be break a huge amount of existing
> code.
>
> Just curious, does OSX take longer because you have to deal with
> multiple possible GDAL versions?
>
> Thanks again,
> Andy
>
> On Mon, Nov 2, 2015 at 1:21 PM, Roger Bivand  wrote:
> > Hi Andy,
> >
> > On Mon, 26 Oct 2015, Andy Teucher wrote:
> >
> >> Hi Roger,
> >>
> >> Sorry I haven't been much help here - it's not for lack of trying,
> >> just lack of skills. I've been poring through OGR_write.cpp and trying
> >> to figure out language and the flow of things... a few thoughts
> >> occurred to me (which are probably obvious to you):
> >>
> >
> > Your example was very helpful, the problem was that when writeOGR() and
> the
> > underlying compiled code was written, the model used was that of the ESRI
> > Shapefile driver, rather than OGC SFS, which was first published in 2004.
> > Most older formats (and their drivers in OGR) use a OGR geometry factory
> > function to organize polygons on reading - but this is not present in
> > shapelib (bundled with the R maptools package), so when designing polygon
> > classes for sp, and writing rgdal::writeOGR() about 10 years ago, this
> was
> > not known, and was not revisited subsequently, until you happily spotted
> the
> > bug.
> >
> > GEOS in rgeos uses OGC SFS definitions (the ones you are using below),
> but
> > the SpatialPolygons/Polygons/Polygon structure was already established by
> > then, and a fix (using comment attributes in Polygons objects) was used,
> > rather than a major change to Polygons objects (in 2010).
> >
> > The revision to rgdal::writeOGR() now branches on the input objects - if
> > they have comments (applied in code, not user-facing), these are trusted,
> > and the MultiPolygon objects constructed containing one or more Polygon
> > objects with one exterior ring and zero or more interior rings on that
> > basis. However, if the input objects do not have comment attributes,
> > OGRGeometryFactory::organizePolygons() (C++) is used to create properly
> > organized MultiPolygon objects for export.
> >
> > This means that the pre-OGC SFS drivers continue to work, but that the
> OGC
> > SFS drivers, like GeoJSON, now also work correctly when the second ring
> in
> > an object is not assumed to be an interior ring (rgdal < 1.1-1). A test
> set
> > for all the possible cases has been introduced.
> >
> > The revision is now in source package form on CRAN:
> >
> > https://cran.r-project.org/web/packages/rgdal/index.html
> >
> > and a Windows binary package will be available later on Tuesday I hope.
> The
> > OSX revision may take a little longer, but we are grateful to have what
> we
> > have anyway.
> >
> > Thanks for a clear and very helpful bug report.
> >
> > Best wishes,
> >
> > Roger
> >
> >
> >> 1) Is this bug because holes are standalone polygons in sp objects,
> >> while they are members of polygon objects in geojson, wkt, etc?
> >>
> >> 2) On a related note, I notice that a single polygon with a hole is
> >> identified as a multipolygon containing a single polygon with a hole,
> >> when using writeOGR and one of the affected drivers.  Practically I
> >> don't think this matters too much, but it's not quite right.
> >>
> >> 3) In lines 116-160 the polygon/multipolygon and line/multiline
> >> determination is made for the entire layer. Can/should this be made
> >> for each feature in the layer, as you could have a layer with a
> >> mixture of single polygons and multipolygons?
> >>
> >> 4) I feel like there ought to be one more layer of looping through the
> >> polyons. e.g.:
> >>
> >> loop over each feature;
> >> determine whether each feature is a single polygon (length 1 or
> >> all rings but one are holes), or a multipolygon:
> >>if a polygon loop over and add ring(s);
> >>if a multipolygon loop over polygons:
> >>for each polygon add ring(s) and somehow assign holes to
> >> parent polygons
> >>
> >> Again, apologies if I'm being totally obvious (or way off base).
> >> Thanks for all your work on these packages, and on this issue.
> >>
> >> Andy
> >>
> >>
> >> On Thu, Oct 22, 2015 at 11:07 PM, Roger Bivand 
> >> wrote:
> >>>
> >>> On Fri, 23 Oct 2015, Andy Teucher wrote:
> >>>
>  On my Mac, I am running rgdal 1.0-7 and GDAL 1.11.3
> 
> 
>  sessionInfo()
> 
> >>>
> >>> OK, thanks.
> >>>
> >>> The problem also affects the SQLite driver (the CRAN binaries do not
> have
> >>> this 

Re: [R-sig-Geo] Exporting Multipolygon geojson with rgdal::writeOGR

[Roger Bivand]

On Mon, 2 Nov 2015, Andy Teucher wrote:


Amazing (and fast) work Roger, and thanks for the great explanation.
I can see why you wouldn't want to redo the way Polygons objects are
structured - I'm sure it would be break a huge amount of existing
code.

Just curious, does OSX take longer because you have to deal with
multiple possible GDAL versions?


No, the builds of the underlying external dependencies for the different 
platforms are handled by the CRAN administrators, Simon Urbanek for OSX, 
and I believe he has to trigger a fresh build manually. Windows (Uwe 
Ligges and Brian Ripley) also has the R Winbuilder service, which is worth 
knowing about, and is more automated.


Once the Windows build train compilers are upgraded, we'll move towards 
getting GDAL 2.0.* on these platforms.


For those using OSX, following the R-sig-mac list is very informative - 
see the thread in October on El Capitan ... of course working around 
things like this take time and attention.


Roger



Thanks again,
Andy

On Mon, Nov 2, 2015 at 1:21 PM, Roger Bivand  wrote:

Hi Andy,

On Mon, 26 Oct 2015, Andy Teucher wrote:


Hi Roger,

Sorry I haven't been much help here - it's not for lack of trying,
just lack of skills. I've been poring through OGR_write.cpp and trying
to figure out language and the flow of things... a few thoughts
occurred to me (which are probably obvious to you):



Your example was very helpful, the problem was that when writeOGR() and the
underlying compiled code was written, the model used was that of the ESRI
Shapefile driver, rather than OGC SFS, which was first published in 2004.
Most older formats (and their drivers in OGR) use a OGR geometry factory
function to organize polygons on reading - but this is not present in
shapelib (bundled with the R maptools package), so when designing polygon
classes for sp, and writing rgdal::writeOGR() about 10 years ago, this was
not known, and was not revisited subsequently, until you happily spotted the
bug.

GEOS in rgeos uses OGC SFS definitions (the ones you are using below), but
the SpatialPolygons/Polygons/Polygon structure was already established by
then, and a fix (using comment attributes in Polygons objects) was used,
rather than a major change to Polygons objects (in 2010).

The revision to rgdal::writeOGR() now branches on the input objects - if
they have comments (applied in code, not user-facing), these are trusted,
and the MultiPolygon objects constructed containing one or more Polygon
objects with one exterior ring and zero or more interior rings on that
basis. However, if the input objects do not have comment attributes,
OGRGeometryFactory::organizePolygons() (C++) is used to create properly
organized MultiPolygon objects for export.

This means that the pre-OGC SFS drivers continue to work, but that the OGC
SFS drivers, like GeoJSON, now also work correctly when the second ring in
an object is not assumed to be an interior ring (rgdal < 1.1-1). A test set
for all the possible cases has been introduced.

The revision is now in source package form on CRAN:

https://cran.r-project.org/web/packages/rgdal/index.html

and a Windows binary package will be available later on Tuesday I hope. The
OSX revision may take a little longer, but we are grateful to have what we
have anyway.

Thanks for a clear and very helpful bug report.

Best wishes,

Roger



1) Is this bug because holes are standalone polygons in sp objects,
while they are members of polygon objects in geojson, wkt, etc?

2) On a related note, I notice that a single polygon with a hole is
identified as a multipolygon containing a single polygon with a hole,
when using writeOGR and one of the affected drivers.  Practically I
don't think this matters too much, but it's not quite right.

3) In lines 116-160 the polygon/multipolygon and line/multiline
determination is made for the entire layer. Can/should this be made
for each feature in the layer, as you could have a layer with a
mixture of single polygons and multipolygons?

4) I feel like there ought to be one more layer of looping through the
polyons. e.g.:

loop over each feature;
determine whether each feature is a single polygon (length 1 or
all rings but one are holes), or a multipolygon:
   if a polygon loop over and add ring(s);
   if a multipolygon loop over polygons:
   for each polygon add ring(s) and somehow assign holes to
parent polygons

Again, apologies if I'm being totally obvious (or way off base).
Thanks for all your work on these packages, and on this issue.

Andy


On Thu, Oct 22, 2015 at 11:07 PM, Roger Bivand 
wrote:


On Fri, 23 Oct 2015, Andy Teucher wrote:


On my Mac, I am running rgdal 1.0-7 and GDAL 1.11.3


sessionInfo()



OK, thanks.

The problem also affects the SQLite driver (the CRAN binaries do not have
this driver). Could someone please check the PostGIS driver - my guess is
that all OGC SFS compliant drivers may be affected.

Could somebody 

Re: [R-sig-Geo] Adding spatial tables to existing SpatiaLite DBs

[Roger Bivand]

On Tue, 3 Nov 2015, Barry Rowlingson wrote:


I can create a SpatiaLite DB file and put a layer in it, but if I try
and add another layer, rgdal fails. Example:

Versions etc:

> require(rgdal)
Loading required package: rgdal
Loading required package: sp
prgdal: version: 1.0-7, (SVN revision 559)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
Path to GDAL shared files: /usr/share/gdal/1.11
Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
Path to PROJ.4 shared files: (autodetected)
Linking to sp version: 1.2-0

Create a simple points data set, write it:

> pts = data.frame(x=runif(10),y=runif(10),z=1:10)
> coordinates(pts)=~x+y
> writeOGR(pts, "tmpfile.db", "pts", driver="SQLite",
dataset_option="SPATIALITE=YES")

Note the use of the dataset_option to make this a proper SpatiaLite,
and not just an SQLite table. The output file is about 4Mb and has a
lot of metadata tables in it. The file loads into QGIS which
recognises it as a SpatiaLite table and I can plot the points.

Now try and create another spatial table (this time, "pts2")  in the
same database file:

> writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite",
dataset_option="SPATIALITE=YES")
Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite",
dataset_option = "SPATIALITE=YES") :
 Creation of output file failed


No reported error in 1.1-1, but the original layer is overwritten. The 
layer in the existing dsn is overwritten. Appending is not something 
writeOGR knows about:


writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite", 

dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
writeOGR(pts, "tmpfile.db", "pts2", driver="SQLite", 

dataset_option="SPATIALITE=YES", layer_option="FORMAT=SPATIALITE")
Error in writeOGR(pts, "tmpfile.db", "pts2", driver = "SQLite", 
dataset_option = "SPATIALITE=YES",  :

  layer exists, use a new layer name

ogrListLayers("tmpfile.db")

[1] "pts2"
attr(,"driver")
[1] "SQLite"
attr(,"nlayers")
[1] 1

For now, system("ogr2ogr ..."), or something from gdalUtils? Contribution 
to writeOGR() to implement appending as in 
http://www.gdal.org/drv_sqlite.html?


Roger



I've tried various other options with no success. Is this possible?
Or do I ditch writeOGR for this and create the table by hand... ick...

Barry

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: roger.biv...@nhh.no

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Re: [R-sig-Geo] Help with as.numeric(rownames(over(SpatialPoints(set1, set2, returnList= TRUE)[[1]]))

[chris english]

Hi Alan;
As you say, runs most of the time. I took the liberty of cleaning out the
>'s, removed the call to plyr as it doesn't seem to be used, and the
rm(list=ls()) since I wasn't playing in a sandbox. How frequently does it
misbehave?

I'm sorry I can't offer anything more constructive than the clean up.  I'm
interested in identifying self-avoiding random walks SAW which means I'll
first have to figure out how to implement one, then figure out how to test
for one. And means getting my head around your  "which_next <-
sample(c("bb","dd","ff","hh"),1)" logic.

Chris

## reproducible example code
 #rm(list = ls())
 library(deldir)
 library(sp)
 #library(plyr)
 side_length = 100
 ## Create random SET of XY coordinates (size = 100x100)
 set.seed(11)
 df = data.frame(matrix(sample(1:100,16,replace=TRUE),nrow=8))
 ## Convert df to SPatialPointsDataFrame
 spdf <- SpatialPointsDataFrame(df,df)
 ## deldir() function creates tesselation (voronoi) plot
 z <- deldir(df,plotit=TRUE,wl='te')
 ## tile.list() creates a list of data for tiles
 zz <- tile.list(deldir(df,plotit=TRUE,wl='te'))
 ## Voronoi Polygons Function
 voronoipolygons = function(layer) {
   require(deldir)
   crds = layer@coords
   z = deldir(crds[,1], crds[,2])
   w = tile.list(z)
   polys = vector(mode='list', length=length(w))
   require(sp)
   for (i in seq(along=polys)) {
 pcrds = cbind(w[[i]]$x, w[[i]]$y)
 pcrds = rbind(pcrds, pcrds[1,])
 polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
   }
   SP = SpatialPolygons(polys)
   voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
 y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
 function(x) slot(x, 'ID'
 }
 ## Generate SpatialPolygonsDataFrame to use as input for over() function
 vpl <- voronoipolygons(spdf)
 ## Random Walk Function generates North, South East or West movement
 ## with transit from across screen (like PacMan, going out one side,
 ## coming back on the other side) to prevent getting stuck in corner
 random_walk <- function(step_quantity, step_length, plot = FALSE){
   require(ggplot2)

   walker <- data.frame(matrix(c(0,0), nrow = step_quantity, ncol = 3,
 byrow = T))
   names(walker)[1]<-paste("x")
   names(walker)[2]<-paste("y")
   names(walker)[3]<-paste("which")

   ## Seed random initial starting point
   walker[1,1:2] <- matrix(sample(1:100,2,replace=TRUE),nrow=1)
   walker[1,3] <- as.numeric(rownames(over(SpatialPoints(walker[1,1:2]
),vpl,returnList=TRUE)[[1]]))

  where_to <- as.numeric()

   for(i in 2:step_quantity){
 where_to <- as.numeric()
 where_to <- walker[i-1,1:2]
 which_next <- sample(c("bb","dd","ff","hh"),1)

 if (which_next == "bb") {
   if (walker[i-1,2] == side_length) {where_to[1,2] <- 0
   } else {where_to[1,2] <- walker[i-1,2]+step_length}
 }

 else if (which_next == "dd") {
   if (walker[i-1,1] == 0 ) {where_to[1,1] <- side_length
   } else {where_to[1,1] <- walker[i-1,1]-step_length}
 }

else if (which_next == "ff") {
   if (walker[i-1,1] == side_length) {where_to[1,1] <- 0
   } else {where_to[1,1] <- walker[i-1,1]+step_length}
 }
 else {
   if (walker[i-1,2] == 0) {where_to[1,2] <- side_length
   } else {where_to[1,2] <- walker[i-1,2]-step_length}
 }

 walker[i,1:2] <- where_to
   }

   walker[i,3] <- as.numeric(rownames(over(SpatialPoints(walker[i,1:2]),
   vpl,returnList= TRUE)[[1]]))


   if(plot){
   require(ggplot2)
   p <- ggplot(walker, aes(x = x, y = y))
   p <- p + geom_path()
   print(p)
   }

   return(walker)
 }
 try(transits <- random_walk(5000,1),silent=F)

On Sun, Nov 1, 2015 at 1:35 PM, Alan Briggs  wrote:

> Hello.
>
> Below is a fully repeatable R-Script that I'm having trouble with.
> Generally, here's what I'm trying to do:
>
> 1) Randomly generate a tile.list()
> 2) Randomly generate a new point
> 3) Identify which polygon in the tile.list the new randomly generated point
> is in
>
> This works fine MOST of the time. However, occasionally, I get an error
> returned:
>
> Error in `[<-.data.frame`(`*tmp*`, list, 3, value = numeric(0)) :
> >   replacement has length zero
>
>
> While troubleshooting, I realized I get numeric(0) returned for certain
> sets of new random points when I run the command
> as.numeric(rownames(over(SpatialPoints(walker[i,1:2]),vpl,returnList=
> TRUE)[[1]])). I thought maybe this was a boundary issue, but the points
> don't lie on the edge, nor are they the centroid.
>
> Any help you can provide would be greatly appreciated!
>
> Regards,
>
> Alan
>
> R-Script Below:
>
> ### Help Question for  r-sig-geo ###
> > rm(list = ls())
> > library(deldir)
> > library(sp)
> > library(plyr)
> > side_length = 100
> > ## Create random SET of XY coordinates (size = 100x100)
> > set.seed(11)
> > df = data.frame(matrix(sample(1:100,16,replace=TRUE),nrow=8))
> > ## Convert df to SPatialPointsDataFrame
> > spdf <- 

Re: [R-sig-Geo] Help with as.numeric(rownames(over(SpatialPoints(set1, set2, returnList= TRUE)[[1]]))

[chris english]

Alan,

Possibly a nonsensical way of testing but I ran

rep(try(transits <- random_walk(5000,1),silent=F), 100)

and it completed without issue.

Chris

On Mon, Nov 2, 2015 at 5:33 AM, chris english  wrote:

> Hi Alan;
> As you say, runs most of the time. I took the liberty of cleaning out the
> >'s, removed the call to plyr as it doesn't seem to be used, and the
> rm(list=ls()) since I wasn't playing in a sandbox. How frequently does it
> misbehave?
>
> I'm sorry I can't offer anything more constructive than the clean up.  I'm
> interested in identifying self-avoiding random walks SAW which means I'll
> first have to figure out how to implement one, then figure out how to test
> for one. And means getting my head around your  "which_next <-
> sample(c("bb","dd","ff","hh"),1)" logic.
>
> Chris
>
> ## reproducible example code
>  #rm(list = ls())
>  library(deldir)
>  library(sp)
>  #library(plyr)
>  side_length = 100
>  ## Create random SET of XY coordinates (size = 100x100)
>  set.seed(11)
>  df = data.frame(matrix(sample(1:100,16,replace=TRUE),nrow=8))
>  ## Convert df to SPatialPointsDataFrame
>  spdf <- SpatialPointsDataFrame(df,df)
>  ## deldir() function creates tesselation (voronoi) plot
>  z <- deldir(df,plotit=TRUE,wl='te')
>  ## tile.list() creates a list of data for tiles
>  zz <- tile.list(deldir(df,plotit=TRUE,wl='te'))
>  ## Voronoi Polygons Function
>  voronoipolygons = function(layer) {
>require(deldir)
>crds = layer@coords
>z = deldir(crds[,1], crds[,2])
>w = tile.list(z)
>polys = vector(mode='list', length=length(w))
>require(sp)
>for (i in seq(along=polys)) {
>  pcrds = cbind(w[[i]]$x, w[[i]]$y)
>  pcrds = rbind(pcrds, pcrds[1,])
>  polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
>}
>SP = SpatialPolygons(polys)
>voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1],
>  y=crds[,2], row.names=sapply(slot(SP, 'polygons'),
>  function(x) slot(x, 'ID'
>  }
>  ## Generate SpatialPolygonsDataFrame to use as input for over() function
>  vpl <- voronoipolygons(spdf)
>  ## Random Walk Function generates North, South East or West movement
>  ## with transit from across screen (like PacMan, going out one side,
>  ## coming back on the other side) to prevent getting stuck in corner
>  random_walk <- function(step_quantity, step_length, plot = FALSE){
>require(ggplot2)
>
>walker <- data.frame(matrix(c(0,0), nrow = step_quantity, ncol = 3,
>  byrow = T))
>names(walker)[1]<-paste("x")
>names(walker)[2]<-paste("y")
>names(walker)[3]<-paste("which")
>
>## Seed random initial starting point
>walker[1,1:2] <- matrix(sample(1:100,2,replace=TRUE),nrow=1)
>walker[1,3] <- as.numeric(rownames(over(SpatialPoints(walker[1,1:2]
> ),vpl,returnList=TRUE)[[1]]))
>
>   where_to <- as.numeric()
>
>for(i in 2:step_quantity){
>  where_to <- as.numeric()
>  where_to <- walker[i-1,1:2]
>  which_next <- sample(c("bb","dd","ff","hh"),1)
>
>  if (which_next == "bb") {
>if (walker[i-1,2] == side_length) {where_to[1,2] <- 0
>} else {where_to[1,2] <- walker[i-1,2]+step_length}
>  }
>
>  else if (which_next == "dd") {
>if (walker[i-1,1] == 0 ) {where_to[1,1] <- side_length
>} else {where_to[1,1] <- walker[i-1,1]-step_length}
>  }
>
> else if (which_next == "ff") {
>if (walker[i-1,1] == side_length) {where_to[1,1] <- 0
>} else {where_to[1,1] <- walker[i-1,1]+step_length}
>  }
>  else {
>if (walker[i-1,2] == 0) {where_to[1,2] <- side_length
>} else {where_to[1,2] <- walker[i-1,2]-step_length}
>  }
>
>  walker[i,1:2] <- where_to
>}
>
>walker[i,3] <- as.numeric(rownames(over(SpatialPoints(walker[i,1:2]),
>vpl,returnList= TRUE)[[1]]))
>
>
>if(plot){
>require(ggplot2)
>p <- ggplot(walker, aes(x = x, y = y))
>p <- p + geom_path()
>print(p)
>}
>
>return(walker)
>  }
>  try(transits <- random_walk(5000,1),silent=F)
>
> On Sun, Nov 1, 2015 at 1:35 PM, Alan Briggs  wrote:
>
>> Hello.
>>
>> Below is a fully repeatable R-Script that I'm having trouble with.
>> Generally, here's what I'm trying to do:
>>
>> 1) Randomly generate a tile.list()
>> 2) Randomly generate a new point
>> 3) Identify which polygon in the tile.list the new randomly generated
>> point
>> is in
>>
>> This works fine MOST of the time. However, occasionally, I get an error
>> returned:
>>
>> Error in `[<-.data.frame`(`*tmp*`, list, 3, value = numeric(0)) :
>> >   replacement has length zero
>>
>>
>> While troubleshooting, I realized I get numeric(0) returned for certain
>> sets of new random points when I run the command
>> as.numeric(rownames(over(SpatialPoints(walker[i,1:2]),vpl,returnList=
>> TRUE)[[1]])). I thought maybe this was a boundary issue, but the points
>> don't lie on the edge,