Hello, I've solved my own problem. I thought I would post my solution in case others find it useful as my task relates to: 1) mosiacing rasters of different extents, 2) working with the PRISM, climateSA, climateNA and climateWNA datasets, 3) replacing values in a raster based on xy coordinates...
To recap, I had: a) ânorthAmâ = one large raster of all of North America (originally from climateNA); b) âdfmâ = a data frame listing just the x y coordinates of (non-NA) cells from a) that correspond to Mexico, as well as the variable variable that I wish to map (e.g. after extracting data from climateSA); columns are x,y, clim1 in that order; c) âusaâ = a raster with the same extent etc. as a) but with values of interest (clim1) for just the USA and NA values everywhere else (e.g. PRISM data). My goal was to produce a final raster with the extent of a) but with values from b) and c). NOTE: I found it easier to work with the data frame of Mexico rather than a raster that can be produced using dataframe2asc (SDMTools). I did the following: # First turn the northAm into a template/empty rasterâ¦will keep all the properties of a) which I want northAm[]<-NA # get the cell IDs for Mexico from the North American raster but based on XY from the Mexico data frame (because I initially derived dfm from the northAm raster, my xy values in dfm directly correspond to those in northAm)  cellID<-cellFromXY(northAm,dfm[,1:2]) # where columns 1 and 2 in dfm are x and y respectively dfm$cellID<-cellID # create a new data frame with just cell IDs and the values from clim1 dfm_new<-dfm[,c(4,3)] # where column 4 is cellID and column 3 is clim1 # now replace values in the template raster with clim1 from Mexico northAm[dfm_new$cellID]<-dfm_new$clim1 # And finally merge this data with the data from the USA final<-merge(northAm,usa) Cheers, Julie ________________________________ From: "[email protected]" <[email protected]> To: [email protected] Sent: Thursday, July 3, 2014 12:00:03 PM Subject: R-sig-Geo Digest, Vol 131, Issue 3 Send R-sig-Geo mailing list submissions to    [email protected] To subscribe or unsubscribe via the World Wide Web, visit    https://stat.ethz.ch/mailman/listinfo/r-sig-geo or, via email, send a message with subject or body 'help' to    [email protected] You can reach the person managing the list at    [email protected] When replying, please edit your Subject line so it is more specific than "Re: Contents of R-sig-Geo digest..." Today's Topics:  1. Help mosiacing rasters OR replacing NA values in one   raster    based on cell coordinates (Julie Lee-Yaw)  2. Re: line color in plot(SpatialPolygonsDataFrame) [was: maps    showing post-Soviet states?] (Spencer Graves)  3. spatial.tools 1.4.8 now on CRAN (Jonathan Greenberg)  4. Create SpatialPolygonsDataFrame indicating the polygons in a    single source SpatialPolygonsDataFrame that overlap with the new    polygons? (nevil amos)  5. Re: Choropleth map (Janka VANSCHOENWINKEL) ---------------------------------------------------------------------- Message: 1 Date: Wed, 2 Jul 2014 05:37:11 -0700 To: "[email protected]" <[email protected]> Subject: [R-sig-Geo] Help mosiacing rasters OR replacing NA values in    one   raster based on cell coordinates Message-ID:    <[email protected]> Content-Type: text/plain Hi Sorry for the complicated example, but I'm really stuck: I am trying to merge two rasters/datasets: one for the continental USA and one for Mexico. The dataset for Mexico was initially derived from a larger raster of North America that was at the same extent, resolution etc. as the US raster...all in an equal area projection (so the Mexico raster simply had NA values for the US and vice versa). For purposes of populating the Mexico layer with appropriate cell values (using a separate program), I needed to convert the Mexico layer to points, reproject the points to WGS84 and feed those points into a separate program using the following: m<-as.data.frame(mexico,xy=TRUE) m_noNA<-na.omit(m) # reproject to WGS84, feed resulting points into external program (climateSA)... My data frame output from the external program includes x and y coordinates for each cell based on the original equal area projection as well as cell values for the variable of interest. I can either leave this output as a data frame OR convert it to an ascii grid using: library(SDMTools) dataframe2asc(myTable) But this asci will have a different extent from the original mexico layer because I had to remove the NA values from the data frame (no way to avoid this). So now I'm trying to figure out how to either: 1) merge/mosaic two rasters with different extents and NO overlap (e.g. the US raster has NA values in Mexico and the Mexico raster [of a different extent] has NA values for the US...). In this case, I don't know how to get merge to produce a raster with cell values for BOTH the US and Mexico (I always end up with NA values in one region or the other). And similarly with mosaic, I don't know what function to use ...there shouldn't be overlap between there rasters and if there is by chance, I would want to populate the resulting raster with the values from the US (not the mean, min, etc.) or 2) For the US raster, replace NA values that overlap with Mexico with the values in myTable based on the x y co-ordinates...I haven't been able to figure out this type of indexing/ look-up yet... [[elided Yahoo spam]] Thanks Julie    [[alternative HTML version deleted]] ------------------------------ Message: 2 Date: Wed, 02 Jul 2014 08:11:12 -0700 From: Spencer Graves <[email protected]> To: Michael Sumner <[email protected]> Cc: "[email protected]" <[email protected]> Subject: Re: [R-sig-Geo] line color in plot(SpatialPolygonsDataFrame)    [was: maps showing post-Soviet states?] Message-ID: <[email protected]> Content-Type: text/plain Hi, Michael:    Thanks for the reply. Your suggestions fixed my plot problem AND my inability to find S4 plot methods. [The latter is now documented with help('methods'), but it wasn't when I last looked for it a few years ago.]    Best Wishes,    Spencer Graves On 7/1/2014 9:11 PM, Michael Sumner wrote: > Here's one way, cast to Spatial*Lines and control the object colours > directly: > > ct1 <- subset(ct, NAME %in% c("Azerbaijan", "Uzbekistan", "Tajikistan", >              "Turkmenistan", "Croatia")) > plot(as(ct1, "SpatialLinesDataFrame"), >      col=seq_len(nrow(ct1)), lwd=3) > > > You can get the source of these S4 methods like this: > > showMethods("plot") > > that tells you the base signature (x="SpatialPolygons", y="missing") > (or SpatialLines), which can be found with > > findMethods("plot")[["SpatialPolygons#missing"]] > > which tells you the source is in the (unexported) function which > ultimately traverses the objects to be plotted with lower level > plotting functions: > > sp:::plot.SpatialPolygons > > Cheers, Mike. > > > > On Wed, Jul 2, 2014 at 1:54 PM, Spencer Graves > <[email protected] <mailto:[email protected]>> wrote: > >  How can one control the color of lines with >  plot(SpatialPolygonsDataFrame)? > > >       Adding lwd=9 to the arguments for plot produces black >  borders so wide most of the red disappears: > > >  map(xlim=c(10, 90), ylim=c(30, 60)) >  library(raster) >  ct <- getData("countries") >  plot(subset(ct, NAME %in% c("Azerbaijan", "Uzbekistan", "Tajikistan", >                "Turkmenistan", "Croatia")), >        add=TRUE, col='red', lwd=9) > > >       A related question is how can I see the plot function >  used here? Methods dispatch does not take this to plot.default nor >  to plot.SpatialPolygonsDataFrame, as it would with an S3 method. > > >       Thanks, >       Spencer > > >  ######################### >     Thanks to Gilles Leduc, Greg Snow, Thomas Adams, Barry >  Rowlingson, and Michael Sumner for their replies. For the >  archives, in case someone else might find this in the archives, I >  will record here two solutions to my problem. > > >     1. Michael's reply was the simplest, especially since it >  worked with "map" adding "add=TRUE" and "col" to the plot: > >  map(xlim=c(10, 90), ylim=c(30, 60)) >  library(raster) >  ct <- getData("countries") >  plot(subset(ct, NAME %in% c("Azerbaijan", "Uzbekistan", "Tajikistan", >                "Turkmenistan", "Croatia")), >        add=TRUE, col='red') > > >     2. Before I saw Michael's reply, I had solved the problem >  with much greater effort: (1) I downloaded a zip file for each >  country from "www.diva-gis.org/datadown >  <http://www.diva-gis.org/datadown>", as suggested by Gilles. (2) >  Unzipping produced folders with names like "..\TWN_adm" and files >  with names beginning "TWN_adm0", "TWN_adm1", and "TWN_adm2" and >  extensions bdf, prj, sbn, sbx, shp, and shx. I could then get >  what I wanted with code like the following: > > >  Taiwan0 <- readShapeSpatial("..\TWN_adm0") >  map(xlim=c(100, 150), ylim=c(10, 30)) >  plot(Taiwan0, add=TRUE, col='blue', lwd=3, fill=FALSE) > > >     Again, thanks for the replies. > > >  On 7/1/2014 4:38 PM, Michael Sumner wrote: > >    An example: > >    library(raster) >    ct <- getData("countries") >    plot(subset(ct, NAME %in% c("Azerbaijan", "Uzbekistan", >    "Tajikistan", >    "Turkmenistan", "Croatia")) > > > > >    On Wed, Jul 2, 2014 at 9:09 AM, Michael Sumner >    <[email protected] <mailto:[email protected]>> wrote: > > >      On Wed, Jul 2, 2014 at 7:56 AM, Thomas Adams >      <[email protected] <mailto:[email protected]>> wrote: > >        Spencer, > >        Maybe I'm not sure what you're asking but here >        http://www.gadm.org/version1 >        are current shapefiles by country. The highlighting >        (selection) canals be >        done in GRASS GIS, QGIS, GMT > > >       raster::getData can read this. See here >      http://www.gadm.org/ and > >      library(raster) >      ?getData > > > > > > >        Tom > > >        On Tue, Jul 1, 2014 at 2:29 PM, Spencer Graves < >        [email protected] >        <mailto:[email protected]>> >        wrote: > >              What would you suggest I do to >create a >          world map highlighting >          Taiwan and post-Soviet states like Azerbaijan, >          Uzbekistan, Tajikistan, >          Turkmenistan, and Croatia? > > >              The region argument in maps::map >won't >          recognize any of these, and >          my literature search identified hundreds of >          packages with some spatial / >          geographical / mapping content. > > >              Thanks, >              Spencer Graves > >  _______________________________________________ >  R-sig-Geo mailing list >  [email protected] <mailto:[email protected]> >  https://stat.ethz.ch/mailman/listinfo/r-sig-geo > > > > > -- > Michael Sumner > Software and Database Engineer > Australian Antarctic Division > Hobart, Australia > e-mail: [email protected] <mailto:[email protected]>    [[alternative HTML version deleted]] ------------------------------ Message: 3 Date: Wed, 2 Jul 2014 17:51:05 -0500 From: Jonathan Greenberg <[email protected]> To: "[email protected]" <[email protected]>,    "[email protected]"    <[email protected]>,    Carlos Ramirez    <[email protected]>,   George Scheer <[email protected]> Subject: [R-sig-Geo] spatial.tools 1.4.8 now on CRAN Message-ID:    <CABG0rfvNJWp9+RJZer4nVTFEwJh+C75hOdYXgqN-g9GiJYof=a...@mail.gmail.com> Content-Type: text/plain; charset=UTF-8 Folks: I released a new version of spatial.tools on CRAN (it should be making its rounds via update.packages() in the next 24 hours). In addition to the requisite bug squashing, here are some the major new features/functions: rasterEngine/focal_hpc: - new parameter: "datatype" -- will now output files using a user-specific datatype (e.g. you can now select integer outputs instead of being forced to use floating point). Check ?dataType in the raster package for a list of possible datatypes. - new parameter: "clearworkers" (TRUE by default) -- once rasterEngine is complete, it will now "flush" any parallel workers registered with foreach (recovering memory). If you have any issues with this [[elided Yahoo spam]] - Useful bug fix: you should get an ENVI header corresponding to the output file before rasterEngine begins its main processing run, so you can monitor the output in semi-real time by simply opening the output file in your favorite GIS program (black stripes are yet-to-be-completed regions). New functions: - raster_to_GLT: Creates a geometry lookup (GLT) file from a raster (a two-band brick where the pixel values for the first band are the column numbers, and the pixel values for the second band are the row numbers of the corresponding pixels in the input Raster* file) - raster_to_IGM: Creates an input geometry (IGM) file from a raster (a two-band brick where the pixel values for the first band are the geographic x-coordinates, and the pixel values for the second band are the geographic y-coordinates of the corresponding pixels in the input Raster* file) [[elided Yahoo spam]] --j -- Jonathan A. Greenberg, PhD Assistant Professor Global Environmental Analysis and Remote Sensing (GEARS) Laboratory Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 259 Computing Applications Building, MC-150 605 East Springfield Avenue Champaign, IL 61820-6371 Phone: 217-300-1924 http://www.geog.illinois.edu/~jgrn/ AIM: jgrn307, MSN: [email protected], Gchat: jgrn307, Skype: jgrn3007 ------------------------------ Message: 4 Date: Thu, 3 Jul 2014 16:17:23 +1000 From: nevil amos <[email protected]> To: [email protected] Subject: [R-sig-Geo] Create SpatialPolygonsDataFrame indicating the    polygons in a single source SpatialPolygonsDataFrame that overlap with    the new polygons? Message-ID:    <CAN9eD7k5Di=bRnVy=NzmJ=vouebrznezkxzafv+pea81un4...@mail.gmail.com> Content-Type: text/plain I am trying to create a polygon dataset that creates a new polygon for each unique combination of location and overlap of polygons in an input polygon dataset, and is attributed with the input polygons underlying each of the derived polygons, and th output does not contain any overlapping polygons. #trivial example #single layer with three partially overlapping polygons library(sp) p1<-Polygons(list(Polygon(cbind(x=c(0,2,2,0,0),y=c(0,0,2,2,0)))),1) p2<-Polygons(list(Polygon(cbind(x=c(1,3,3,1,1),y=c(1,1,3,3,1)))),2) p3<-Polygons(list(Polygon(cbind(x=c(1.5,3,3,1.5,1.5),y=c(1.5,1.5,3,3,1.5)))),3) mypolys<-list(p1,p2,p3) mydata<-data.frame(Name=c("p1","p2","p3"),row.names=c(1,2,3)) SDF<-SpatialPolygonsDataFrame(SpatialPolygons(mypolys),mydata) par(mfrow=c(3,4)) plot(SDF) #the desired result is a SpatialPolygonsDataFrame with a single separate polygon for each of the six areas deifined by the interection of the boundaries of the three input polygons text(.5,.5,1) text(1.25,1.25,2) text(1.75,1.75,3) text(2.5,2.5,4) text(1.25,2.5,5) text(2.5,1.25,6) # such that the data frame lists the input polygons (p1,p2,p3) that overlap at the location of each of the output polygons (1:6) and the output has no overlap between the polygons # the results would be the six unique polygons shown in the plot and so the attributes would look like: data.frame(polys=c("p1","p1 p2","p1 p2 p3","p2 p3","p2","p2")) # or alternatively data.frame(cbind(c("p1","p1","p1","p2","p2","p2"),c(NA,"p2","p2","p3",NA,NA),c(NA,NA,"p3",NA,NA,NA))) # or alternatively data.frame(cbind(p1=c(1,1,1,0,0,0),p2=c(0,1,1,1,1,1),p3=c(0,0,1,1,0,0))) #This last plot shows the "apparent polygons" I want in separate colors ( because you cannot see the overlapping portions) Outpoly<-union(SDF,SDF) plot(OutPoly,col=1:9) plot(SpatialPolygons(list(Polygons(list(Polygon(cbind(x=c(2,2,3,3,2),y=c(1,1.5,1.5,1,1)))),1))),col=12,add=T) text(.5,.5,1) text(1.25,1.25,2) text(1.75,1.75,3) text(2.5,2.5,4) text(1.25,2.5,5) text(2.5,1.25,6) #The nearest I could get was raster::union however this requires two ( and only 2) SpatialDataFrames as input , what I need is union with a single input union(SDF,SDF) results in 9 polygons not 6, these are all simple rectangles of overlapping pairs (union(union(SDF,SDF)SDF)) results in 27 polygons - all again just the simple rectangular overlaps and original polygons union(SDF,SDF)@data    [[alternative HTML version deleted]] ------------------------------ Message: 5 Date: Thu, 3 Jul 2014 10:35:58 +0200 From: Janka VANSCHOENWINKEL <[email protected]> To: James Rooney <[email protected]>, [email protected] Cc: "[email protected]" <[email protected]> Subject: Re: [R-sig-Geo] Choropleth map Message-ID:    <CAHymut+20D_QqSyJs=zzeq0sojhabyqqttcm3n3akt63+gc...@mail.gmail.com> Content-Type: text/plain Hi James and Roger, You both pointed out that maybe I should not use the merge function to join data to an sp object. I found the explanation for this on the internet ( http://gis.stackexchange.com/) and thought it would be interesting to tell you: The merge function resorts the data during the operation which breaks the internal relationship in the sp object. Therefore to merge a dataframe to the @data slot of an sp object, using match is indeed better and is done in this way: shape = name of shapefile Otherdata = data that I want to match IDS = identifier I want to merge on shape@data = data.frame(shape@data, OtherData[match(shape@data$IDS, OtherData$IDS),]) Furthermore, I have not yet entirely found the code to make the chloropleth [[elided Yahoo spam]] Thanks again for you advice, Janka 2014-06-25 14:40 GMT+02:00 James Rooney <[email protected]>: > Dear Janka, > > I can't answer your ggplot questions as I don't really use it much however > I did spot something in your code that will cause problems. > > You used a "merge" to combine data between a dataframe and an SPDF: > > > ValueMapDf <- merge(nuts, Value, by.x="NUTS_ID", by.y="NUTS_ID") > > I've forgotten the reasoning, but this will not give you the merge results > you expect. You might want to do some manual checking on your merge results. > James > > ________________________________________ > From: [email protected] [[email protected]] > On Behalf Of Janka VANSCHOENWINKEL [[email protected]] > Sent: 24 June 2014 18:35 > To: [email protected] > Subject: [R-sig-Geo] Choropleth map > > Dear list members, > > > I am trying to make a basic choropleth map for all nuts3 regions in Europe > ( = a map where > > regions are colored by the value of a certain variable ââ¬â think for > instance > about a populations > > density map). > > > I am struggling with this for over a week by now and I always get the same > error message: > > Error: Aesthetics must either be length one, or the same length as > > the dataProblems:ValueMapDf$NUTS_ID > > > > I checked the book ââ¬Åggplot2â⬠from Hadley Weckham and several > internet > sources. They > > document the different steps very properly, but still I donââ¬â¢t know > what I > am doing wrong. > > [[elided Yahoo spam]] > > *Below you can find my code. I hope I gave sufficient information, > otherwise, just let me know!* > > > library(maptools) > > library(ggplot2) > > library(ggmap) > > > # read administrative boundaries (I call the SpatialPolygonsDataFrame: > "nuts" since it contains all nuts3 regions of Europe) ) > > nuts <- readShapePoly(fn="NUTS_RG_60M_2006") > > names(nuts) > > [1] "OBJECTID" "NUTS_ID"  "STAT_LEVL_" "AREA"   "LEN" > "Shape_Leng" "Shape_Area" > > > # Then I have a dataset with over 60000 observations (several observations > per nuts3 region). I take the average of the different > > observations per nuts3 region and save it as "Value" > > # The dataset is called Alldata and MEt is the variable I would like to map > over the different nuts3 regions > > > Value <- aggregate(Alldata$MEt,list(Alldata$nuts3),mean) > > colnames(Value) <- c("NUTS_ID","ValueMEt") > > > head(Value) > > NUTS_ID ValueMEt > > 1 AT111 0.04856170 > > 2 AT112 0.05448523 > > 3 AT113 0.04563415 > > 4 AT121 0.02164469 > > 5 AT122 0.02664611 > > 6 AT123 0.03163034 > > > > # merge map and data: I call the result "ValueMapDf" > > ValueMapDf <- merge(nuts, Value, by.x="NUTS_ID", by.y="NUTS_ID") > > ValueMapDf <- ValueMapDf[order(ValueMapDf$NUTS_ID),] > > > > # then I dropped some regions of the SpatialPolygonsDataFrame, for which I > don't have observations (I did this manually since I didn't > > know if there was a way to do it automatically) > > ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID > 441,] > > ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 444,] > > ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 447,] > > ââ¬Â¦. > > ValueMapDf<-ValueMapDf[ValueMapDf$OBJECTID != 1927,] > > > [[elided Yahoo spam]] > However, as you can see below, here something goes wrong. > > # ggplot mapping > > # data layer > > m0 <- ggplot(data=ValueMapDf) > > # fill with Value data > > m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=ValueMapDf$NUTS_ID, > fill="ValueMapDf$ValueMEt")) + coord_equal() > > # add map with only borders (to have visible borders) > > m2 <- m1 + geom_path(aes(x=long, y=lat, group=ValueMapDf$NUTS_ID), > color='black') > > m2 > > Error: Aesthetics must either be length one, or the same length as the > dataProblems:ValueMapDf$NUTS_ID > > > > I also tried: > > qplot(long,lat,data=ValueMapDf, group=ValueMapDf$NUTS_ID, > fill=ValueMapDf$ValueMEt, geom="polygon") > > Regions defined for each Polygons > > Error: Aesthetics must either be length one, or the same length as the > dataProblems:ValueMapDf$ValueMEt, ValueMapDf$NUTS_ID > > > > > > > length(ValueMapDf$OBJECTID) > > [1] 1122 > > > length(ValueMapDf) > > [1] 1122 > > > length(ValueMapDf$NUTS_ID) > > [1] 1122 > > > length(ValueMapDf$ValueMEt) > > [1] 1122 > > > > *When I write:* > > m0 <- ggplot(data=ValueMapDf) > > # fill with Value data > > m1 <- m0 + geom_polygon(aes(x=long, y=lat, group=NUTS_ID, > fill="ValueMapDf$ValueMEt")) + coord_equal() > > # add map with only borders (to have visible borders) > > m2 <- m1 + geom_path(aes(x=long, y=lat, group=NUTS_ID), color='black') > > m2 > > It get: Error in eval(expr, envir, enclos) : object 'NUTS_ID' not found > > > [[elided Yahoo spam]] > > > Janka > >    [[alternative HTML version deleted]] > >    [[alternative HTML version deleted]] ------------------------------ _______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo End of R-sig-Geo Digest, Vol 131, Issue 3 ***************************************** [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list [email protected] https://stat.ethz.ch/mailman/listinfo/r-sig-geo
