






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: r.hijmans@gmail.com
To: sglish@hotmail.com
CC: r-sig-geo@r-project.org


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 <sglish@hotmail.com> 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

R-sig-Geo@r-project.org

https://stat.ethz.ch/mailman/listinfo/r-sig-geo




 		 	   		  
	[[alternative HTML version deleted]]

