issue in using projectRaster in raster library
Responses below:
On Thu, Apr 10, 2014 at 4:49 PM, ping yang <pingyang.whu at gmail.com> wrote:
Hi Jonathan, I am using GDAL 1.9.2 and -wo NUM_THREADS=ALL_CPUS didn't work. Actually I have several thousands file to process, so I am thinking of the foreach, I changed the code as follows:
GDAL 1.9.2 is pretty old (October 2012) -- I'd recommend upgrading to 1.10.1 -- I think they added the parallel support around v. 1.10.0. OSGeo4W is my install of choice for Windows, but you can see the various versions at: http://trac.osgeo.org/gdal/wiki/DownloadingGdalBinaries Also, you need to use the -multi parameter as well, otherwise it will ignore the -wo NUM_THREADS=ALL_CPUS.
ProjectANDSlope<- function(zipfile)
{
require(raster)
require(rgdal)
require(gdalUtils)
require(foreach)
zipfile <-filelist[i]
#First Unzip the file
unzip(zipfile,overwrite = TRUE,exdir = "../UNZIPPED")
longitude <- paste("-",substr(zipfile,5,7),sep="")
lon <- as.numeric(longitude)
zoneNo <- (floor((lon + 180)/6) %% 60) + 1
utm <-paste("+proj=utm +zone=",zoneNo," +datum=NAD83 +ellps=GRS80
+towgs84=0,0,0 +units=m +no_defs",sep="")
outputfile <- paste("../PROJECTED/",substr(zipfile,1,7),".tif",sep='')
inputfile <- paste("../UNZIPPED/grd",substr(zipfile,1,7),"_13",sep="")
#Project Raster to UTM
gdalwarp(srcfile=inputfile, dstfile=outputfile,
t_srs=utm,output_Raster=TRUE,overwrite=TRUE,multi=TRUE,
tr=c(10,10),r="bilinear", verbose=TRUE)
#Generate Slope raster
slopefile <- paste("../SLOPE/",substr(zipfile,1,7),".tif",sep='')
gdaldem(mode="slope",input_dem=outputfile,
output=slopefile,compute_edges=TRUE,output_Raster=TRUE,verbose=TRUE)
#delete the intermedia files
unlink(inputfile,recursive=TRUE,force=TRUE)
}
and use different methods to run it, First I use:
filelist <-dir(".",pattern="n*.zip") #all the file list
system.time(
for(i in 21:30)
{
ProjectANDSlope(filelist[i])
}
)
The time used
user system elapsed
25.58 36.31 810.08
Then I used foreach (do and dopar)
system.time(
foreach(i=21:30) %dopar% ProjectANDSlope(filelist[i])
)
The time used:
user system elapsed
24.86 35.50 757.93
not much improved, then I realized that I didn't generate the parallel
backend, I added the following code:
require(parallel)
Sys.getenv('NUMBER_OF_PROCESSORS') "8"
cl <-makeCluster(4)
But I get errors:
Error in socketConnection("localhost", port = port, server = TRUE, blocking
= TRUE, :
cannot open the connection
In addition: Warning message:
In socketConnection("localhost", port = port, server = TRUE, blocking =
TRUE, :
port 11216 cannot be opened
This looks like a potential firewall problem -- but you also need to
register it with foreach -- this is a good question for r-sig-hpc, by
the way, since that is where the foreach people tend to hang out. try
something like:
cl <- makeCluster(spec=4,type="PSOCK")
registerDoParallel(cl)
Keep in mind, if you don't see your CPUs being pegged it is because
you are (likely) I/O limited.
Incidentally, the "cleaner" way to use packages with foreach is using
the .packages= parameter rather than using require():
foreach(...,.packages=c("raster","rgdal","gdalUtils")) # You don't
need to require foreach -- it is autoloaded
I don't know what happened(I remember I used the makeCluster before but now it is not working). By the way, can you help to specify how to compare the two methods ( parallel processing each individual image vs. parallel processing multiple images at the same time)?
The idea is for a lot of files, you can approach it in two ways:
1) Sequentially loop through each file, but parallelize the single
file processing (e.g. using -multi -wo NUM_THREADS=ALL_CPUS in GDAL)
-- when you see all the processors light up, this is processing ONE
file.
2) Parallel loop through your files -- each CPU ("worker"), then, is
processing a SINGLE file, but you are processing multiple files at the
same time. This would be more like what you are trying to accomplish
above (using foreach).
Slope calculations are a form of focal-window analysis which can be
CPU limited, so I think #2 is the right way to go. If you are trying
to cut down the time, you should look for bottlenecks -- are you
reading and writing to the same drive? Are you doing it over a
network? Things like that can really slow the process down.
--j
Look forward to hearing from you. Thanks, Ping On Tue, Apr 8, 2014 at 5:42 PM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Ping: Re: raster:::projectRaster I'm not sure how much you can speed that up (Robert?) -- I think it is native R code which is generally a lot slower than compiled code like the GDAL utilities. It has to use "chunking", but for reasonable chunk sizes this is not likely to cause any noticeable overhead. Re: gdalUtils vs. cygwin -- you got basically the same execution time -- 21.15 vs. 20.7 seconds is not worth worrying about (~ a 2% increase using gdal from R)-- gdalUtils has the (small) overhead from R, but it is actually using your installed gdalwarp in cygwin anyway, it is just an R wrapper for it if you are more comfortable using R functions vs. command line. The first time you launch gdalwarp in an R session it takes a few seconds to find your local installation of GDAL -- the second time it runs, there should be almost no overhead. Re: parallelization -- if you have GDAL 1.10.1 or later (perhaps a slightly earlier version, not sure when it was implemented), you can use -multi -wo NUM_THREADS=ALL_CPUS to parallelize the operation -- in gdalUtils, this would be: gdalwarp(...,multi=TRUE,wo="NUM_THREADS=ALL_CPUS") # See: http://lists.osgeo.org/pipermail/gdal-dev/2012-June/033084.html If you are trying to warp a LOT of images, you can also use e.g. foreach() to warp each image using a different thread. You'd have to experiment with which one (loop through each image one at a time, parallel processing each individual image vs. parallel processing multiple images at the same time) is more efficient. --j On Tue, Apr 8, 2014 at 4:50 PM, ping yang <pingyang.whu at gmail.com> wrote:
Hi Jonathan,
I am following up with the previous thread, Now I successfully installed
the
gdalUtils from the source and ran successfully the task (both projection
and
slope calculation) using the following script:
rasterOptions(chunksize = 3e+07, maxmemory = 1e+08) #suggestions
according
to Jonanthan, it made the raster::projectRaster work but for much longer
time such as 40 minutes!!!
#*----Function for detect the UTM Zone-----#
long2UTM <- function(long) {
(floor((long + 180)/6) %% 60) + 1
}
filename <- "grdn41w083_13"
longitude <- paste("-",substr(filename,9,10),sep="") #this could be
change
according to the directory string
lon <- as.numeric(longitude)
utm <-paste("+proj=utm +zone=",long2UTM(lon)," +datum=NAD83
+ellps=GRS80
+towgs84=0,0,0 +units=m +no_defs",sep="")
outputfile <-
paste("projected/",substr(filename,4,10),"_raster.tif",sep='')
gdalwarp(srcfile=filename, dstfile=outputfile,
t_srs=utm,output_Raster=TRUE,overwrite=TRUE,multi=TRUE,
tr=c(10,10),r="bilinear", verbose=TRUE)
slp <- gdaldem(mode="slope",input_dem=outputfile,
output="n41w083_rgdalslp.tif",output_Raster=TRUE,verbose=TRUE)
Now it took me 21.15 seconds to finish the projection(while using the
gdalWarp in cygwin took about 20.7 seconds)
I am wondering if there is way to speed up this process, How to parallel
this process(parallelized warping)?
Big Thanks,
Ping
On Mon, Apr 7, 2014 at 11:37 PM, Jonathan Greenberg <jgrn at illinois.edu>
wrote:
For gdalUtils, you don't need to upgrade to 3.2 -- the latest CRAN version only requires R 2.14 or later. It should find your cygwin version, but let me know if it doesn't. The key with the rasterOptions you used is that, previously, you were forcing the ENTIRE RASTER to load into memory with the max memory and chunk size you used -- this is really not a good idea for 99.9% of raster calculations. I would try two things: 1) reduce the max memory back to its default settings (1e+08 cells on my machine -- about 750mb of RAM usage if each cell is 64 bits), and 2) adjust the chunk size to maybe 250mb (3e+07 cells) If you leave your settings like that, particularly on Windows, you'll end up with out-of-memory crashes fairly often (since I'm guessing you don't have 8TB of RAM on a Windows box). --j On Mon, Apr 7, 2014 at 4:47 PM, ping yang <pingyang.whu at gmail.com> wrote:
Hi Alex and Jonathan, The reason I used the rasterOptions is because when previously I did the reproject for 30 meter DEM, it worked after I raised the chunksize(otherwise it ran out of memory), so I added here just to avoid that issue(probably it not work for this case) I am going to install gdalUtils package, however, I need to upgrade my R from 3.0.1 to 3.2 (I am afraid to do it since some other work is running). However, I have an update for this, i used the gdalwarp in my Cygwin(actually installed by the QGIS) and run it successfully: gdalwarp -t_srs '+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0' -r bilinear -multi -tr 10 10 -of GTiff grdn41w083_13/hdr.adf n41w083.tif it took 11 second to finish. I am still looking for the slop calculation based on the projected tif file, I will appreciated if anyone could give me some suggestions if have this experience. Thanks a lot. Ping On Mon, Apr 7, 2014 at 3:42 PM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
Ping:
Re: maxmemory -- you really shouldn't adjust that, definitely not to
what you set it to -- you've set it to (by my quick calculations) to
be around 7500 GB of needed memory (this is assuming each cell takes
64 bits of space) -- this is probably why you are crashing.
Chunksize
also is usually set correctly, raster processes generally have a
quickly diminishing returns as you increase the chunk size.
You might want to try gdalwarp in my gdalUtils package -- it will be
faster to reproject, and you can even enable parallelized warping.
install.packages("gdalUtils",
repos="http://R-Forge.R-project.org",type="source") # To get the
latest, otherwise just use CRAN
--j
On Mon, Apr 7, 2014 at 3:30 PM, Alex Zvoleff
<azvoleff at conservation.org>
wrote:
On Mon, Apr 7, 2014 at 4:05 PM, ping yang <pingyang.whu at gmail.com> wrote:
Dear r-sig-geo, I plan to reproject from lonlat to UTM for a bunch of 10 meter NED files (large) using a R script(I did it successfully for the 30 meter DEM), I used the following code to run: rasterOptions(chunksize = 1e+10, maxmemory = 1e+12) #need to allow bigger chunksize for the processing!!!
Why do you find that you "need to allow bigger chunksize"? I rarely need to adjust the chunksize, even for large images.
#*----Function for detect the UTM Zone-----#
long2UTM <- function(long) {
(floor((long + 180)/6) %% 60) + 1
}
require(raster)
require(rgdal)
filename <- "grdn41w083_13"
longitude <- paste("-",substr(filename,9,10),sep="")
lon <- as.numeric(longitude)
utm <-paste("+proj=utm +zone=",long2UTM(lon)," +datum=NAD83
+ellps=GRS80
+towgs84=0,0,0 +units=m +no_defs",sep="")
r <- raster(filename)
outputfile <-
paste("projected/",substr(filename,4,10),".tif",sep='')
print(paste(outputfile))
pr <- projectRaster(r, crs=utm, res=10.0, method='bilinear')
However, the projectRaster takes forever to finish and it almost
crashed my
computer.
The file is kind of big (file size is 455MB),The file
information
is
class : RasterLayer
dimensions : 10812, 10812, 116899344 (nrow, ncol, ncell)
resolution : 9.259259e-05, 9.259259e-05 (x, y)
extent : -83.00056, -81.99944, 39.99944, 41.00056 (xmin,
xmax,
ymin,
ymax)
coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+no_defs
data source : c:\temporary\10M\grdn41w083_13
names : grdn41w083_13
values : 212.254, 460.0799 (min, max)
My R information:
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
other attached packages:
[1] rgdal_0.8-10 raster_2.1-49 sp_1.0-11
loaded via a namespace (and not attached):
[1] grid_3.0.1 lattice_0.20-27 tools_3.0.1
However, this process only took 10 seconds using the existing
tool
in
ArcGIS toolbox. Has someone experienced this before? Any comments
will
be
welcome.
Are you using "Project" or "Define Projection" in ArcGIS? "Define projection" just overwrites the existing projection rather than reprojecting the data, which is much more calculation intensive.
Thanks,
Ping
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Alex Zvoleff Postdoctoral Associate Tropical Ecology Assessment and Monitoring (TEAM) Network Conservation International 2011 Crystal Dr. Suite 500, Arlington, Virginia 22202, USA Tel: +1-703-341-2749, Fax: +1-703-979-0953, Skype: azvoleff http://www.teamnetwork.org | http://www.conservation.org On Mon, Apr 7, 2014 at 4:05 PM, ping yang <pingyang.whu at gmail.com> wrote:
Dear r-sig-geo,
I plan to reproject from lonlat to UTM for a bunch of 10 meter
NED
files
(large) using a R script(I did it successfully for the 30 meter
DEM), I
used the following code to run:
rasterOptions(chunksize = 1e+10, maxmemory = 1e+12) #need to
allow
bigger
chunksize for the processing!!!
#*----Function for detect the UTM Zone-----#
long2UTM <- function(long) {
(floor((long + 180)/6) %% 60) + 1
}
require(raster)
require(rgdal)
filename <- "grdn41w083_13"
longitude <- paste("-",substr(filename,9,10),sep="")
lon <- as.numeric(longitude)
utm <-paste("+proj=utm +zone=",long2UTM(lon)," +datum=NAD83
+ellps=GRS80
+towgs84=0,0,0 +units=m +no_defs",sep="")
r <- raster(filename)
outputfile <-
paste("projected/",substr(filename,4,10),".tif",sep='')
print(paste(outputfile))
pr <- projectRaster(r, crs=utm, res=10.0, method='bilinear')
However, the projectRaster takes forever to finish and it almost
crashed my
computer.
The file is kind of big (file size is 455MB),The file
information
is
class : RasterLayer
dimensions : 10812, 10812, 116899344 (nrow, ncol, ncell)
resolution : 9.259259e-05, 9.259259e-05 (x, y)
extent : -83.00056, -81.99944, 39.99944, 41.00056 (xmin,
xmax,
ymin,
ymax)
coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0
+no_defs
data source : c:\temporary\10M\grdn41w083_13
names : grdn41w083_13
values : 212.254, 460.0799 (min, max)
My R information:
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
attached base packages:
[1] stats graphics grDevices utils datasets methods
base
other attached packages:
[1] rgdal_0.8-10 raster_2.1-49 sp_1.0-11
loaded via a namespace (and not attached):
[1] grid_3.0.1 lattice_0.20-27 tools_3.0.1
However, this process only took 10 seconds using the existing
tool
in
ArcGIS toolbox. Has someone experienced this before? Any comments
will
be
welcome.
Thanks,
Ping
[[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Alex Zvoleff Postdoctoral Associate Tropical Ecology Assessment and Monitoring (TEAM) Network Conservation International 2011 Crystal Dr. Suite 500, Arlington, Virginia 22202, USA Tel: +1-703-341-2749, Fax: +1-703-979-0953, Skype: azvoleff http://www.teamnetwork.org | http://www.conservation.org [[alternative HTML version deleted]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- 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: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
-- 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: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
-- 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: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
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: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007