Home
May 2, 2018 - Mike McMahon (mike.mcmahon@dfo-mpo.gc.ca)
This package contains functions that are either required by and/or useful to other packages from http://github.com/Maritimes. Most of the functions originated in other packages but were moved here to facilitate their use by other packages.
At this time, data with privacy considerations must be aggregated such that each polygon has a minimum of 5 unique values for sensitive fields like Licenses, License Holders, and Vessels. This function takes a dataframe and shapefile and for each polygon in the shapefile calculates:
- Aggregate values for a number of (user-specified) fields , and
- How many unique values exist in each polygon for each of a number of sensitive fields.
The following two shapefiles can be generated:
- A shapefile looks like the original submitted shapefile, but where every polygon is populated with:
- The original shapefile fields;
- The MEAN, COUNT and SUM of all of the fields identified as
agg.fields; - One or more columns showing the number of unique values for each of the fields identified as
sens.fields; - A column ("TOTUNIQUE"), indicating the number of unique value for the least represented
sens.field. For example, if a polygon has 10 unique values for LICENSE, but only 4 for VRN, this column would display "4". The purpose is to be as conservative as possible an ensure that privacy is maintained. - A column ("CAN_SHOW") indicating whether or not more detailed, gridded data can be shown from that polygon. Values of YES and NO show that data was present, while blank fields mean there were no records in the polygon to assess.
- A shapefile of 2 min grid cells, where only the data that lay within polygons where "CAN_SHOW = Yes" is gridded. No gridded data at all is available for areas where "CAN_SHOW = No".
assess_privacy(df=data, agg.fields = c("TOTWGT", "TOTNO"),agg.poly.field = "NAFO_BEST", sens.fields = c("MISSION", "SETNO"), create.shps = T)Since no shapefile was supplied as agg.poly.shp, NAFOSubunits gets used instead.
| Example Image | Description |
|---|---|
 |
As an example, imagine the data on the left had sensitivity issues (it's actually non-sensitive RV survey data for redfish). Now pretend that "MISSION" and "SETNO" are sensitive, personally-identifiable fields (like "LICENCE_NUMBER" and "VRN"). We want to share the data, but we have to make sure that we have enough different values of "MISSION" and "SETNO" in any given area so that no one can figure out which data is attributable to which "MISSION" or "SETNO". |
 |
If we run the code shown above this table, we get two shapefiles. The original data points are still there, albeit faded. Also shown is "screened_areas*.shp", and it is symbolized by "CAN_SHOW". Open areas have "CAN_SHOW=YES" Orange, hashed areas have "CAN_SHOW=NO" Greyed areas had no data. Polygon labels show the different NAFO areas, and in brackets, the number of unique, private records within each (i.e. "TOTUNIQUE"). |
 |
On the left is a snippet of the attributes of the "screened_areas*.shp" shown above. In it, you can see: The MEAN, COUNT and SUM values for each of the agg.fields; the number of unique values for each of the sens.fields;"TOTUNIQUE", showing the minimum number of all of the sens.fields;and"CAN_SHOW", showing whether or not data from within this polygon can be aggregated and displayed. |
 |
The second shapefile created is "2MinGrid*.shp". This is the original point data gridded to 2 min squares for those areas where data is allowed to be shown (i.e. "CAN_SHOW = Y" in screened_areas*.shp). In th eexample on the left, the grid cells are coloured by the field "TOTNO_SUM", one of the fields calculated from an agg.field. Note that no gridded data can be seen in the areas that were shown "CAN_SHOW=NO" (i.e. the orange, hashed areas). |
 |
Zooming in on the "2MinGrid*.shp" shows how grid cells are not even generated for those polygon areas where "CAN_SHOW = N". |
 |
The table on the left shows the fields that are available for each grid cell - again, note that the fields available are calculated from the fields identified in agg.field. |
This function takes a dataframe and converts in to a SpatialPointsDataframe.
The sp package cannot plot a data.frame (with valid coordinates) by default:
> str(data)
'data.frame': 1435 obs. of 5 variables:
$ MISSION : chr "ATC1970176" "ATC1970176" "ATC1970176" "ATC1970176" ...
$ SETNO : int 10 13 15 17 1 23 26 28 2 32 ...
$ LATITUDE : num 44.7 44.8 44.2 44.3 44.3 ...
$ LONGITUDE: num -60.9 -60.1 -59.9 -60.5 -62.2 ...
$ TOTNO : int 1642 1107 41 560 947 7 3463 102 369 246 ...
> sp::plot(data)
Error in plot.window(...) : need finite 'xlim' values
...but if we run the same data against this function, it gets converted from a data.frame to a SpatialPointsDataFrame, and can be plotted.
> data_sp = df_to_sp(data)
> str(data_sp)
Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
...
>sp::plot(data_sp)
This function identifies points that aren't in the northern hemisphere (i.e. LAT between 0 and 90) and aren't in the western hemisphere (i.e. LON between 0 and 180). By default, it returns the "good" points, but the "return.bad" flag allows it to return only the "bad" points. Southern hemisphere points are not technically invalid, but for Maritimes pruposes, they should never come up.
In the example below, we see that by default, running df_qc_spatial on our data gives 1435 records. If we manually make 2 of the coordinates bad, and run it again, those 2 records are omitted. If we want to see what the records were, we can set return.bad=TRUE and they are revealed.
#check how many rows we find by default
> nrow(df_qc_spatial(data))
[1] 1435
# manually mess up 2 records
> data[1,"LATITUDE"]<-91
> data[2,"LONGITUDE"]<-NA
> nrow(df_qc_spatial(data))
[1] 1433
# we find 2 less coordinates
>sp::plot(df_to_sp(data))
Error in .local(obj, ...) : NA values in coordinates
# and we can't plot them!
# if we do our spatial qc prior to plotting them, we can!
>data_clean = df_qc_spatial(data)
>sp::plot(df_to_sp(data_clean))
>
# let's see which records are bad
> df_qc_spatial(data, return.bad = T)
MISSION SETNO LATITUDE LONGITUDE TOTNO
1 ATC1970176 10 91.00000 -60.88333 1642
2 ATC1970176 13 44.81667 NA 1107
> This function takes a dataframe and removes any columns that are entirely populated with NAs
> head(data)
MISSION SETNO LATITUDE LONGITUDE TOTNO GARBAGE
1 ATC1970176 10 91.00000 -60.88333 1642 NA
2 ATC1970176 13 44.81667 NA 1107 NA
3 ATC1970176 15 44.16667 -59.85000 41 NA
4 ATC1970176 17 44.33333 -60.51667 560 NA
5 ATC1970176 1 44.31667 -62.20000 947 NA
6 ATC1970176 23 43.41667 -60.56667 7 NA
>data_clean = drop_NACols(data)
> head(data_clean)
MISSION SETNO LATITUDE LONGITUDE TOTNO
1 ATC1970176 10 91.00000 -60.88333 1642
2 ATC1970176 13 44.81667 NA 1107
3 ATC1970176 15 44.16667 -59.85000 41
4 ATC1970176 17 44.33333 -60.51667 560
5 ATC1970176 1 44.31667 -62.20000 947
6 ATC1970176 23 43.41667 -60.56667 7
This function takes a df with coordinate fields in decimal degrees, and overlays it a shapefile (agg.poly.shp). It adds a column to the df indicating which polygon within the shapefile each point falls within. If no polygon is provided, the df will be assessed against NAFO subdivisions.
> #default uses NAFO Subdivisions
> head(identifyArea(df=data))
MISSION SETNO LATITUDE LONGITUDE TOTNO NAFO_BEST
3 ATC1970176 15 44.16667 -59.85000 41 <NA>
4 ATC1970176 17 44.33333 -60.51667 560 4WE
5 ATC1970176 1 44.31667 -62.20000 947 4WK
6 ATC1970176 23 43.41667 -60.56667 7 4WG
7 ATC1970176 26 43.53333 -59.83333 3463 4WG
8 ATC1970176 28 43.65000 -59.30000 102 4WG
>
> # but we can use custom shapefiles and fields too
> head(identifyArea(df=data, agg.poly.shp = "/mnt/R_PED/Shared/Spatial/Science/Strata/ped_groundfish/MaritimesRegionEcosystemAssessmentStrata(2014-).shp",agg.poly.field = "StrataID"))
MISSION SETNO LATITUDE LONGITUDE TOTNO StrataID
3 ATC1970176 15 44.16667 -59.85000 41 457
4 ATC1970176 17 44.33333 -60.51667 560 457
5 ATC1970176 1 44.31667 -62.20000 947 462
6 ATC1970176 23 43.41667 -60.56667 7 454
7 ATC1970176 26 43.53333 -59.83333 3463 497
8 ATC1970176 28 43.65000 -59.30000 102 453This function takes a df of the isdb data (including fields FISHSET_ID and LAT1:LAT4 and LONG1:LONG4) and makes a spatialLinesDataFrame which can be plotted.If selected, it can also QC them. Each line will have the following columns in the resulting data frame:
- FISHSET_ID This uniquely identifies a set
- QCPOS This field identifies potential issues with the line - including NA positions, positions that are "0", incorrect hemisphere, or impossible coordinates
- QCTIME This field identifies potential issues with the timing of the vertices making up the line. For example, P2 should occur after P1, and 2 vertices cannot occur simultaneously. Sets less than 5 min or greater than 24 hr are also noted
- LEN_KM This field shows the calculated distance of the resultant line in kms
- N_VALID_VERT This field shows how many vertices appear correct after NAs, 0s and othe problematic values have been dropped
#to simply plot the tracks, you can do this
plot(make_isdb_tracks(isdb.df=ISSETPROFILE_WIDE))
# to generate a qc report on the tracks, you can do this
> sets = make_isdb_tracks(isdb.df=ISSETPROFILE_WIDE,do.qc = T, return.choice = "notlines")
#I cheated and found an area where the QC shows some variety, and changed the FISHSET_ID to nonsense
FISHSET_ID QCTIME QCPOS LEN_KM N_VALID_VERT
2219 999996971 No Time between pts Zero_Len Line 0.0000000 NA
2220 999996972 No Time between pts 0.3767921 2
2221 999996973 No Time between pts Zero_Len Line 0.0000000 NA
2222 999996974 No Time between pts Zero_Len Line 0.0000000 NA
2223 999996975 No Time between pts Zero_Len Line 0.0000000 NA
2224 666679489 Zero_Len Line 0.0000000 NA
2225 666679490 Zero_Len Line 0.0000000 NA
2226 666679491 Zero_Len Line 0.0000000 NA
2227 666679492 Zero_Len Line 0.0000000 NA
2228 666679493 Zero_Len Line 0.0000000 NA
2229 666679494 Zero_Len Line 0.0000000 NA
2230 666679495 Zero_Len Line 0.0000000 NA
2231 666679496 Zero_Len Line 0.0000000 NA
2232 666679497 Zero_Len Line 0.0000000 NA
2233 666679398 Zero_Len Line 0.0000000 NA
2234 666682028 0.4586141 2
This function facilitates creating a connection to Oracle, and allows connection via RODBC or ROracle, depending on the value of usepkg. Credentials can be passed directly to the function, but if they're left blank, the function will prompt you.
Note that the connection object itself isn't what's returned, but a list of 3 objects:
| usepkg chosen | test$usepkg | test$channel | test$thecmd |
|---|---|---|---|
test=make_oracle_cxn(usepkg = 'rodbc') |
'rodbc | "RODBC" | sqlQuery |
test=make_oracle_cxn(usepkg = 'roracle') |
'roracle' | "OraConnection" | dbGetQuery |
> test=make_oracle_cxn(usepkg = 'rodbc')
Oracle Username: username
[1] "username"
Oracle Password: mypassword
[1] "mypassword"
Oracle DSN (e.g. PTRAN): PTRAN
[1] "PTRAN"
Successfully connected to Oracle via RODBC
> class(test$channel)
[1] "RODBC"
#lets test it
> test$thecmd(test$channel,"SELECT * FROM dual")
DUMMY
1 X#send the credentials
> test2=make_oracle_cxn(usepkg = 'roracle',fn.oracle.username = "username",fn.oracle.password = "password",fn.oracle.dsn = "PTRAN")
Successfully connected to Oracle via ROracle
> class(test2$channel)
[1] "OraConnection"
attr(,"package")
[1] "ROracle"
#lets test it
> test$thecmd(test$channel,"SELECT * FROM dual")
DUMMY
1 XThis function takes a dataframe with coordinates in decimal degrees and a track identifier, and creates line segments for each distinct identifier. If specified, it will plot the results in R and/or create shapefiles.
While making segments, you can also generate a shapefile of all of the positions/vertices.
-
points = "all"will create a point shapefile that includes all of the positions. -
points = "none"will not generate a point shapefile. -
points = "orphans"will only create points for those "tracks" that only have a single position
> segments = makeSegments(data, objField = "MISSION", seqField = "SETNO",lat.field = "LATITUDE", lon.field = "LONGITUDE")
Shapefiles are limited in the naming conventions of the fields - they can only be 10 characters, and must be devoid of special characters that are allowed in R. This function attempts to ensure that the shortened names that result from calls to writeOGR() can still be understood and are not simply truncated.
> names(data)
[1] "MISSION" "SETNO" "LATITUDE" "LONGITUDE" "NUMBER.OF_FISH%CAught"
> names(prepare_shape_fields(data))
[1] "MISSION" "SETNO" "LATITUDE" "LONGITUDE" "NUMBER_OF_"Converts a vector of characters into something that an "IN" statement can use. By default, the vector values are surround by apostrophes, as though they were holding character values. Setting apos=FALSE, removes these, making the result appropriate for numeric values.
> vec = c(1:10)
> vec
[1] 1 2 3 4 5 6 7 8 9 10
> paste0("AND data IN (",SQL_in(vec),")")
[1] "AND data IN ('1','2','3','4','5','6','7','8','9','10')"
> paste0("AND data IN (",SQL_in(vec,apos = FALSE),")")
[1] "AND data IN (1,2,3,4,5,6,7,8,9,10)"Calculates the standard error of a vector
> vec = c(1:10)
> st_err(vec)
[1] 0.9574271This function extracts VMS data for a given timespan and area. A time buffer can be added returns other points that are not within the area of interest, but give context to the path.
Interestingly, for the example extraction, almost 100 records had bad positions, so wrapping the VMS in a df_qc_spatial() allows us to plot them.
vms = VMSGetRecs(dateStart = '2010-06-15 00:00:00',dateEnd = '2010-06-15 12:00:00',vrnList = NULL )
> nrow(vms)
# makeSegments() is ideal for plotting this kind of data (or turning it into shapefiles)
> vmsSegs = makeSegments(df_qc_spatial(vms), objField = "VR_NUMBER", seqField = "POSITION_UTC_DATE", lat.field = "LATITUDE", lon.field="LONGITUDE")
To get fancier, you can even ask it to only return VMS data that lies within a certain spatial area, and in such cases, it's useful to add an hrBuffer, which adds data from the surrounding hours. This gives some context to what the vessel was doing around the time it was crossing the area of interest. If your area shapefile has multiple polygons, the value for shp.field will be tacked onto your returned data so you know which polygon a particular vessel was in.
In the example below, we'll extract data that existed within the 5Z9 strata on a particular day with a variety of hrBuffer values.
vms5Z9 = VMSGetRecs(dateStart = '2010-06-15 00:00:00',dateEnd = '2010-06-17 00:00:00',shp = "/home/mike/sf_Documents/5Z9.shp",shp.field = "StrataID",hrBuffer = 2)
vms5Z9Segs = makeSegments(vms5Z9, objField = "VR_NUMBER", seqField = "POSITION_UTC_DATE", lat.field = "LATITUDE", lon.field="LONGITUDE")| No hour buffer (hrBuffer = 0) | No hour buffer (hrBuffer = 2) | No hour buffer (hrBuffer = 10) |
|---|---|---|
 |
 |
 |
Additionally, a vector or VRNs can be provided to limit the results to certain vessels. In the interest of privacy, I'm not putting up an example of that, but the if you wanted the tracks of a particular fleet as it related to particular polygons, you might do:
test = VMSGetRecs(dateStart = '2010-06-15 00:00:00',dateEnd = '2010-06-17 00:00:00',shp = "C:/specialAreas.shp",shp.field = "FULL_DESCR",hrBuffer = 2, vrnList =c(11111,11112,11113))
head(test)
VR_NUMBER LATITUDE LONGITUDE POSITION_UTC_DATE SPEED_KNOTS UPDATE_DATE FULL_DESCR SEGMID
11111 42.340000 -63.13150 2010-06-15 00:07:00 NA 2012-03-23 16:47:18 <NA> 1525361819_11111_1
11111 42.340000 -63.13150 2010-06-15 00:07:00 NA 2012-03-23 16:47:18 Special Area 1 1525361819_11111_1
11112 41.350000 -63.13833 2010-06-15 00:17:36 NA 2012-03-23 16:47:21 <NA> 1525361819_11112_1
11112 41.350000 -63.13833 2010-06-15 00:17:36 NA 2012-03-23 16:47:21 Special Area 1 1525361819_11112_1
11112 41.340000 -65.13583 2010-06-15 00:47:53 NA 2012-03-23 16:47:24 Special Area 1 1525361819_11112_1
11113 42.920000 -65.13950 2010-06-15 01:00:00 NA 2012-03-23 16:48:07 <NA> 1525361819_11113_1
# the VR_NUMBERS and positions have all been changed to silly valuesThis is a SpatialPolygonsDataFrame of the NAFO divisions. It is derived from the file that was available on the VDC. The primary fields of interest are NAFO_1, NAFO_2, NAFO_3 and NAFO_BEST. NAFO_BEST is always populated with the most detailed NAFO subunit possible, while the _1, _2 and _3 in the names of the other columns essentially refer to how many characters exist in each field following the integer.
> tail(NAFOSubunits@data)
OBJECTID Shape_Leng Shape_Area NAFO_1 NAFO_2 NAFO_3 NAFO_BEST
0 1 33.508870 5.1751418 6A <NA> <NA> 6A
1 2 36.842889 12.6304564 6B <NA> <NA> 6B
109 192 10.905534 1.0828271 4W 4WD <NA> 4WD
110 193 3.999996 0.8055399 4W 4WE <NA> 4WE
140 249 6.907884 2.5544606 3P 3PS 3PSH 3PSH
141 290 7.354711 1.7175268 4V 4VS 4VSB 4VSBShow code for plotting full dataset.
sp::plot(NAFOSubunits)
Show code for plotting detailed extent (with labels)
# change the extent
NAFOSubunits@bbox[1,]<-c(-63,-58)
NAFOSubunits@bbox[2,]<-c(45,48)
# plot it
sp::plot(NAFOSubunits)
# label it by NAFO_BEST
library(maptools)
invisible(text(maptools::getSpPPolygonsLabptSlots(NAFOSubunits), labels=as.character(NAFOSubunits$NAFO_BEST), cex=1.0))
This is a SpatialPolygonsDataFrame of 2min grid cells for the Maritimes region. This grid is often used when aggregating data.
sp::plot(grid2Min)
A smaller extent shows the detail.
# change the extent
grid2Min@bbox[1,]<-c(-63.80,-63.60)
grid2Min@bbox[2,]<-c(44.40,44.70)
# plot it
sp::plot(grid2Min)
# add the nafo subunits for context
sp::plot(NAFOSubunits, add=T, border="blue") #add the NAFO areas for context and relative size
It should be noted that in the 2 min grid squares were intentionally constructed so that they did *NOT### overlap lines of longitude and latitude. Since data is often reported to the nearest degree or minute, the offset grid ensures that these common values fall into the appropriate cell, and is not subject to "edge effects" that sometime occur when spatial data lands on a line. Below, crosshairs indicate the intersections of the minutes of longitude and latitude.
