using this loop (from FlowingData example) on a subset of cleaned up IRS data
for 40 states where gcIntermediateis provided sufficient data to execute, i.e.
lines are formed and plotted:
## get rid of cases that aren't running all the way throughmovein_0910 <-
unique(inflow0910$st_abbv)movein_0910 <-
movein_0910[!grepl("TX$|IL$|MN$|CO$|NM$|AZ$|CA$|WA$|PR$|AK$|NV$|UT$",
movein_0910)]
for (i in 1:length(movein_0910)) {
map("world", col="#f2f2f2", fill=TRUE, bg="white", lwd=0.05, xlim=xlim,
ylim=ylim)m_insub <- inflow0910[inflow0910$st_abbv == movein_0910[1,],]
##options "NY",] or movein_0910[1],] m_insub <-
m_insub[order(m_insub$return_num),]max_return_num <- max(m_insub$return_num)
for (j in 1:(length)(m_insub$fips_in_orig)) { message(j) ##for(k
in 1:length) (inter)) { ##message(k) movein_orig <-
county_centroid[county_centroid$fips == (m_insub[j,]$fips_in_orig),]
movein_dest <- county_centroid[county_centroid$fips ==
(m_insub[j,]$fips_in_dest),]
inter <- gcIntermediate(c(movein_orig[1,]$long, movein_orig[1,]$lat),
c(movein_dest[1,]$long, movein_dest[1,]$lat),sp=TRUE, n=100,
addStartEnd=TRUE) lines_in <<- inter colindex <-
round((m_insub[j,]$return_num / (max_return_num/10000)) * length(colors))#per
Paul Butler Facebook visualize lines(inter, col=colors[colindex],
lwd=0.6) } ##}
dev.off() }
There is a similar loop for outflow data with similar results
Same loop, case of Texas:
4647Error in movein_dest[1, ] : subscript out of boundsIn addition: Warning
messages:1: In min(x) : no non-missing arguments to min; returning Inf2: In
max(x) : no non-missing arguments to max; returning -Inf3: In min(x) : no
non-missing arguments to min; returning Inf4: In max(x) : no non-missing
arguments to max; returning -Inf
> traceback()5: movein_dest[1, ]4: movein_dest[1, ]3: inherits(p,
> "SpatialPoints")2: .pointsToMatrix(p2)1: gcIntermediate(c(movein_orig[1,
> ]$long, movein_orig[1, ]$lat), c(movein_dest[1, ]$long, movein_dest[1,
> ]$lat), sp = TRUE, n = 100, addStartEnd = TRUE) at #8>
> str(movein_orig)Formal class 'SpatialPointsDataFrame' [package "sp"] with 5
> slots ..@ data :'data.frame': 1 obs. of 7 variables: {snip] ..@
> coords : num [1, 1:2] -96.8 32.8 .. ..- attr(*, "dimnames")=List of 2
> .. .. ..$ : chr "2807" .. .. ..$ : chr [1:2] "coords.x1" "coords.x2" ..@
> bbox : num [1:2, 1:2] -96.8 32.8 -96.8 32.8 .. ..- attr(*,
> "dimnames")=List of 2 .. .. ..$ : chr [1:2] "coords.x1" "coords.x2" .. ..
> ..$ : chr [1:2] "min" "max" [snip]butstr(moveout_dest)Formal class
> 'SpatialPointsDataFrame' [package "sp"] with 5 slots ..@ data
> :'data.frame': 19 obs. of 7 variables: [snip] ..@ coords : num [1:19,
> 1:2] -178 -176 -176 -176 -176 ... .. ..- attr(*, "dimnames")=List of 2 ..
> .. ..$ : chr [1:19] "675" "676" "677" "678" ... .. .. ..$ : chr [1:2]
> "coords.x1" "coords.x2" ..@ bbox : num [1:2, 1:2] -178.3 21.4 -157.8
> 28.4 .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr [1:2] "coords.x1"
> "coords.x2" .. .. ..$ : chr !
[1:2] "min" "max" [snip]in cases that work @data is 1 obs. for 7 variables for
both moveout_orig and moveout_destand @coords $ char "xxx" is the same for both.
1) I would really like to clear up what is happening in my failed cases
2) probably easily overlooked in the above loop the line: line_in <<- interthat
I hope to pass to raster and gdistance to net the lines andand plot the
resultant, perhaps on a net mortality chloropleth...but my efforts this far a
very tentative for this leg.
script, data, shapefiles available at
https://github.com/chris-english/R_GIS
Thank you for your kind attention and advice greatly appreciated
Chris> sessionInfo()R version 2.15.1 Patched (2012-08-15 r60278)Platform:
i386-w64-mingw32/i386 (32-bit)
locale:[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252 [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages: [1] splines tcltk tools stats graphics
grDevices utils datasets [9] methods base
other attached packages: [1] gdistance_1.1-3 Matrix_1.0-6
igraph0_0.5.5-2 raster_2.0-08 [5] maptools_0.8-16
lattice_0.20-6 foreign_0.8-50 rgdal_0.7-12 [9]
Hmisc_3.9-3 survival_2.36-14 sqldf_0.4-6.4
RSQLite.extfuns_0.0.1[13] RSQLite_0.11.1 chron_2.3-42
gsubfn_0.6-4 proto_0.3-9.2 [17] DBI_0.2-5
debug_1.2.40 mvbutils_2.5.101 stringr_0.6.1 [21]
geosphere_1.2-28 sp_0.9-99 maps_2.2-6
loaded via a namespace (and not attached):[1] cluster_1.14.2 grid_2.15.1
plyr_1.7.1Date: Sun, 19 Aug 2012 22:30:45 -0700
Subject: Re: [R-sig-Geo] Great Circles US IRS Migration mapping
From: [email protected]
To: [email protected]
CC: [email protected]
Chris,
I am guessing that you are looking for something like this:
in_start <- !is.na(county$m_0910_in_orig)
m_0910_in_orig <- county[in_start, ]
in_end <- !is.na(county$m_0910_in_dest)
m_0910_in_dest <- county[in_end, ]
To draw great circles, you can create these use the geosphere package (function
gcIntermediate)
Robert
On Sun, Aug 12, 2012 at 9:38 PM, Chris English <[email protected]> wrote:
Hi,
I am trying to map US migration from IRS files using great circles. The IRS
provides two files (either state or county level)
for origin and destination in and out of the state or county. The code thus far
gives me a spatialpolygonsdataframe:
> str(county)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
..@ data :'data.frame': 3258 obs. of 15 variables:
.. ..$ area : num [1:3258] 0.00303 0.13356 0.00267 0.74387 0.93214
...
.. ..$ perimeter : num [1:3258] 0.257 2.05 0.247 5.424 6.036 ...
.. ..$ co99_d00_ : num [1:3258] 1088 1277 157 158 159 ...
.. ..$ co99_d00_i : num [1:3258] 1087 1276 156 157 158 ...
.. ..$ state : Factor w/ 49 levels "01","04","05",..: 38 37 46 25 46
26 22 49 46 46 ...
.. ..$ county : Factor w/ 301 levels "001","003","005",..: 6 52 33 48
5 65 43 25 41 41 ...
.. ..$ name : Factor w/ 1799 levels "Abbeville","Acadia",..: 1704
1094 1486 1387 303 1040 899 1646 1730 1730 ...
.. ..$ lsad : Factor w/ 5 levels "06","08","09",..: 1 1 1 1 1 1 1 1
1 1 ...
.. ..$ lsad_trans : Factor w/ 3 levels "city","County",..: 2 2 2 2 2 2 2
2 2 2 ...
.. ..$ fips : chr [1:3258] "44009" "42091" "53057" "30085" ...
.. ..$ m_0910_in_dest : int [1:3258] NA NA NA NA NA 58798 NA NA NA NA ...
.. ..$ m_0910_in_orig : int [1:3258] NA NA NA NA NA 58803 NA NA NA NA ...
.. ..$ m_0910_out_dest: int [1:3258] 3472 3054 3209 58014 3269 59154 51486
3505 2117 2117 ...
.. ..$ m_0910_out_orig: int [1:3258] 84440 83024 106089 58018 104703 59149
51481 110502 106622 106622 ...
.. ..$ labpt : num [1:3258, 1:2] -71.6 -75.4 -122.6 -105 -120.6 ...
.. ..- attr(*, "data_types")= chr [1:9] "N" "N" "N" "N" ...
..@ polygons :List of 3258
I used the following to extract the long/lat of the label point @labpt to use
as startpoints or endpoints for the great circle lines
depending on whether they were dest(ination) or orig(in) and dest or orig was
not NA.
## reach down to the @labpt level of county using ggplot2 and gpclib to
## get coordinates of county label points to use as start and endpoints
## for make lines function
require(ggplot2)
require(gpclib)
## in any case, using ggplot2 capabilities
county$labpt = coordinates(county)
county$labpt gives me the centroid of all my polygons and I think these follow
the order 1:3258
> county$labpt[1:1,1:2]
[1] -71.57849 41.17788
> county$fips[1]
[1] "44009"
so county$fips[1] "44009" should have its centroid at [1] -71.57849 41.17788
I'm trying to figure out how I am going to make lines between say
county$in_start[6] as startpoint
> county$in_start <- !is.na(county$m_0910_in_orig)
> county$in_start
[1] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
and county$in_end[6] as endpoint
> county$in_end <- !is.na(county$m_0910_in_dest)
> county$in_end
[1] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
I feel like I have most of the pieces but have lost track somewhere. Any
suggestions
greatly appreciated.
Thanks,
Chris
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
