Skip to content

NA x and y values in the population #80

@LHMarshall

Description

@LHMarshall

Reproducible error

# load simulation package
library(dsims)
library(sf)

#Load shapefile Of the region
shapefile.name <-st_read(dsn="DUStrataSeptember.shp")

#Make sure the projection is in equal area if not dssd become unstable.
sf::st_crs(shapefile.name)
# Define WGS 84 / North Pole LAEA Canada (EPSG: 3573)
proj4string<-"+proj=laea +lat_0=90 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
# Project the study area on to a flat plane
projected.strata <- sf::st_transform(shapefile.name, crs = proj4string)

#Make Region
region <- make.region(region.name = "Study.Area",
                      shape = projected.strata)
# There are 15 strata you had only provided 14 names
region <- make.region(region.name = "Study.Area",
                      strata.name = c("A","B","C","D","E","F","G","H","I","J","K","L","M","N", "O"),
                      shape = projected.strata)
##Maybe the mistake is the strata.name? Dsims just overwrite my "Name" column in my shapefile. 
plot(region)

#Creating a Coverage Grid
sf_use_s2(FALSE)
library(lwgeom)
cover<-make.coverage(region, n.grid.points = 1000)
plot(region, cover)

#Design transect
default.design <- make.design(region = region,
                              transect.type = "line",
                              design = "systematic",
                              spacing = c(13500,5500,10000,14000,8500,8500,4000,9500,14000,14000,8000,12000,15000,15000,8000),
                              design.angle = c(05,35,75,75,45,165,0,45,75,135,45,155,10,155,45),
                              edge.protocol = "minus",
                              truncation = 1000, 
                              coverage.grid = cover)
transects <- generate.transects(default.design)
plot(region, transects)

density <- make.density(region = region, 
                        x.space = 5000, 
                        constant = c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1))

#Area of high/low density are based on point density analysis ran on groups 
# Add some areas of higher 
density <- add.hotspot(object = density,
                       centre = c(-442000, -2330000),
                       sigma = 30000,
                       amplitude = 10)
# Add some areas of higher 
density <- add.hotspot(object = density,
                       centre = c(-360000, -2320000),
                       sigma = 10000,
                       amplitude = 10)
# Add some areas of higher 
density <- add.hotspot(object = density,
                       centre = c(-300000, -2300000),
                       sigma = 10000,
                       amplitude = 10)

# Plot the density surface
plot(density,region)


# Create the population description, with a population size N = 3,825 and need to be listed by strata; fixed.N = True the population
#is generated from the values of N otherwise it is generated from the density grid.
pop.desc <- make.population.description(region = region, 
                                        density = density,
                                        N = c(1,821,111,400,470,89,487,5,0,65,130,236,0,0,10),
                                        fixed.N = TRUE)

#Detectability
# Make a simple half normal detection function with a scale parameter of 225
detect.hn <- make.detectability(key.function = "hn",
                                scale.param = 350, 
                                truncation = 1000)
# We can now visualise these detection functions
plot(detect.hn, pop.desc)

#Analysis
ddf.analyses <- make.ds.analysis(dfmodel = list(~1),
                                 key = c("hn", "hr"),
                                 criteria = "AIC",
                                 truncation = 1000)

#Simulation
sim <- make.simulation(reps = 2, 
                       design = default.design,
                       population.description = pop.desc,
                       detectability = detect.hn,
                       ds.analysis = ddf.analyses)

set.seed(138)
survey <- run.survey(sim)

# Running in debug with breakpoint in line 70 of calc.perp.dists.R
pop[1:5,]
   individual         x        y Region.Label scale.param
NA          1        NA       NA            A         350
1           2 -426366.3 -2293717            B         350
2           3 -449451.1 -2291243            B         350
3           4 -472783.7 -2300455            B         350
4           5 -518231.1 -2283110            B         350

I note that stratum A has a population size of 1

Digging into the functions I find that the first attempt to generate the population results in all points (2xN) being located out with the stratum but the checks miss this as a data frame with NAs is created.

image

Metadata

Metadata

Assignees

Labels

No labels
No labels

Type

No type

Projects

No projects

Milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions