I've got a NetLogo model in which each animal occupies a "territory", wherein all patches that belong to the animal are the same color as the animal. I then use the R extension in NetLogo to create a minimum convex polygon (MCP) surrounding each territory and export those polygons as a shapefile. I then import the GIS file back into NetLogo using the GIS extension. My rationale for this is related to this question I posted before (NetLogo - applying values to patches within polygons). That is, I can use gis:intersecting
to give a variable to those patches that fall within the GIS-imported polygons. The process of creating the MCP, exporting, and importing is done at each time step because the territory updates each tick. All works well except that the imported shapefile does not align perfectly with the original polygons. The attached image shows this, where the blue outlines are from the imported shapefile. . I've tried gis:set-world-envelope (list min-pxcor max-pxcor min-pycor max-pycor)
but to no avail. Does anybody know if I'm doing something wrong or if it is error inherent to exporting and then importing a shapefile for which no projection exists? Any help here would be great because having them align would solve several issues, including the previous post. The whole code is pretty long so I've attached some snippets below that are most relevant. Thanks!
extensions [r gis ]
breed [females female]
globals
[
hr-dataset
]
females-own
[
Name
X
Y
]
patches-own
[
is-hr?
]
to setup
clear-all
r:clear
...
ask n-of 5 patches
[
sprout-females 1
[
...
set X (list pxcor)
set Y (list pycor)
]
]
reset-ticks
end
to go
...
expand
calc-homerange
tick
end
to expand
repeat 10
[
ask females
[
move-to target
set X lput pxcor X
set Y lput pycor Y
]
]
end
to calc-homerange
r:eval "library(adehabitatHR)"
r:eval "library(sp)"
r:eval "library(rgdal)"
; create an empty data.frame"
r:eval "turtles <- data.frame()"
; merge the Name, X- and Y-lists of all females to one big data.frame
ask females
[
(r:putdataframe "turtle" "X" X "Y" Y)
r:eval (word "turtle <- data.frame(turtle, Name = '" Name "')")
r:eval "turtles <- rbind(turtles, turtle)"
]
; create SpatialPointsDataFrame
r:eval "spdf <- SpatialPointsDataFrame(turtles[1:2], turtles[3])"
r:eval "homerange <- mcp(spdf, percent = 100)"
r:eval "writeOGR(homerange, '.', 'homerange-rgdal', overwrite_layer = TRUE, driver = 'ESRI Shapefile')"
mark-homeranges
end
to mark-homeranges
clear-drawing
...
set hr-dataset gis:load-dataset "C:/Program Files (x86)/NetLogo 5.0.4/homerange-rgdal.shp"
gis:set-world-envelope (list min-pxcor max-pxcor min-pycor max-pycor)
gis:set-drawing-color blue
gis:draw hr-dataset 2
ask patches gis:intersecting hr-dataset
[
set is-hr? true
]
end
I think Seth is correct that it's an off-by-one error in mapping patch coordinates to world coordinates. It may be related to the bug that was mentioned here. His proposed solution of using
(list min-pxcor - 0.5 max-pxcor + 0.5 min-pycor - 0.5 max-pycor + 0.5)
should work. If not, send me the model and I'll take a look if I have time.Nice screen shot, thanks for providing that, it makes the problem easy to grasp. It looks like the farther from the origin, the worse the discrepancy is, like the error is radiating out from the origin.
I know the GIS extension documentation recommends
gis:set-world-envelope (list min-pxcor max-pxcor min-pycor max-pycor)
, but I have to wonder if it's actually correct.min/max-pxcor/pycor
are the minimum and maximum patch coordinates, not the minimum and maximum turtle coordinates. So for example ifmax-pxcor
is 10, then a turtle's x coordinate can be as high as 10.499999[...], and ifmin-pxcor
is -10, a turtle's x coordinate can be as low -10.5. Thusworld-width
is 21, not 20.Maybe try changing it to
gis:set-world-envelope (list min-pxcor - 0.5 max-pxcor + 0.5 min-pycor - 0.5 max-pycor + 0.5)
and see if that fixes it. Or if it doesn't fix it exactly, try flipping the signs or try 1 instead of 0.5. It sure looks from the screen shot like the problem is an off-by-one or off-by-0.5 error.