[R-sig-Geo] Writing to Large Raster File
Matteo Mattiuzzi
matteo.mattiuzzi at boku.ac.at
Tue Apr 30 18:55:11 CEST 2013
Hi Tom
In fact there is a difference:
see ?dataType why I do not use 'INT4U'. Probably the reason is that writeStart (see ?writeRaster) uses FLT4S as default and adding values in blocks the function doesn't know that FLT4S is not enough where else with out[] <- 1:ncell(out) it knows! I think that the solution.
> ?dataType
> nColVal<-1900
> outFileName <- 'tmp.grd'
>
> #Set up raster
> out <- raster(nrow=20000,ncol=nColVal)#modisRasterAll
>
> ##Get number of colums and blocksize
> outCol <- ncol(out)
> bs <- blockSize(out)
>
> ##Write to file
> out <- writeStart(out,outFileName,overwrite=TRUE,datatype="FLT8S")
> for(i in 1:bs$n){
+ out <-
+ writeValues(out,cellFromRowCol(out,bs$row[i],1) :
+ cellFromRowCol(out,bs$row[i] + bs$nrow[i] -1 , ncol(out)) ,bs$row[i])
+ }
> out <- writeStop(out)
>
> out2 <- raster(out)
> out2[]<- 1:ncell(out2)
>
> res <- out-out2
tail(res)
1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894
19991 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19992 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19993 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19994 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19995 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19996 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19997 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19998 0 0 0 0 0 0 0 0 0 0 0 0 0 0
19999 0 0 0 0 0 0 0 0 0 0 0 0 0 0
20000 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1895 1896 1897 1898 1899 1900
19991 0 0 0 0 0 0
19992 0 0 0 0 0 0
19993 0 0 0 0 0 0
19994 0 0 0 0 0 0
19995 0 0 0 0 0 0
19996 0 0 0 0 0 0
19997 0 0 0 0 0 0
19998 0 0 0 0 0 0
19999 0 0 0 0 0 0
20000 0 0 0 0 0 0
> tail(out)
1881 1882 1883 1884 1885 1886 1887 1888
19991 37982881 37982882 37982883 37982884 37982885 37982886 37982887 37982888
19992 37984781 37984782 37984783 37984784 37984785 37984786 37984787 37984788
19993 37986681 37986682 37986683 37986684 37986685 37986686 37986687 37986688
19994 37988581 37988582 37988583 37988584 37988585 37988586 37988587 37988588
19995 37990481 37990482 37990483 37990484 37990485 37990486 37990487 37990488
19996 37992381 37992382 37992383 37992384 37992385 37992386 37992387 37992388
19997 37994281 37994282 37994283 37994284 37994285 37994286 37994287 37994288
19998 37996181 37996182 37996183 37996184 37996185 37996186 37996187 37996188
19999 37998081 37998082 37998083 37998084 37998085 37998086 37998087 37998088
20000 37999981 37999982 37999983 37999984 37999985 37999986 37999987 37999988
1889 1890 1891 1892 1893 1894 1895 1896
19991 37982889 37982890 37982891 37982892 37982893 37982894 37982895 37982896
19992 37984789 37984790 37984791 37984792 37984793 37984794 37984795 37984796
19993 37986689 37986690 37986691 37986692 37986693 37986694 37986695 37986696
19994 37988589 37988590 37988591 37988592 37988593 37988594 37988595 37988596
19995 37990489 37990490 37990491 37990492 37990493 37990494 37990495 37990496
19996 37992389 37992390 37992391 37992392 37992393 37992394 37992395 37992396
19997 37994289 37994290 37994291 37994292 37994293 37994294 37994295 37994296
19998 37996189 37996190 37996191 37996192 37996193 37996194 37996195 37996196
19999 37998089 37998090 37998091 37998092 37998093 37998094 37998095 37998096
20000 37999989 37999990 37999991 37999992 37999993 37999994 37999995 37999996
1897 1898 1899 1900
19991 37982897 37982898 37982899 37982900
19992 37984797 37984798 37984799 37984800
19993 37986697 37986698 37986699 37986700
19994 37988597 37988598 37988599 37988600
19995 37990497 37990498 37990499 37990500
19996 37992397 37992398 37992399 37992400
19997 37994297 37994298 37994299 37994300
19998 37996197 37996198 37996199 37996200
19999 37998097 37998098 37998099 37998100
20000 37999997 37999998 37999999 38000000
>>> Tom Philippi 04/30/13 6:27 PM >>>
Ted & Matteo--
Is this perhaps a dataType problem? I don't recall the specifics of the
significand of FLT4S but it's roughly 8 digits, so 38000000 may exceed it.
You might try changing the dataType to INT4U to see if the dataType is the
problem.
Tom 2
On Tue, Apr 30, 2013 at 7:22 AM, Matteo Mattiuzzi <
matteo.mattiuzzi at boku.ac.at> wrote:
> Ted you are right!
>
>
> I have tried a few things and I can not explain the behaviour.
> @Robert may the example below helps a little further.
> I have also tried to write into a tif and an img file but it doesn't
> change. Also after reloading the file out <- raster("tmp.grd") the problem
> persists. Starting at row 8831 the values are strange...
> [I have also changed the value generation part, as it is little bit
> clearer to me for reading: cellFromRowCol(out,bs$row[i],1) :
> cellFromRowCol(out,bs$row[i] + bs$nrow[i] -1 , ncol(out)) ]
>
>
> R version 2.15.3 (2013-03-01) -- "Security Blanket"
> Copyright (C) 2013 The R Foundation for Statistical Computing
> ISBN 3-900051-07-0
> Platform: x86_64-pc-linux-gnu (64-bit)
>
>
> R is free software and comes with ABSOLUTELY NO WARRANTY.
> You are welcome to redistribute it under certain conditions.
> Type 'license()' or 'licence()' for distribution details.
>
>
> Natural language support but running in an English locale
>
>
> R is a collaborative project with many contributors.
> Type 'contributors()' for more information and
> 'citation()' on how to cite R or R packages in publications.
>
>
> Type 'demo()' for some demos, 'help()' for on-line help, or
> 'help.start()' for an HTML browser interface to help.
> Type 'q()' to quit R.
>
>
> > library(raster)
> Loading required package: sp
> > sessionInfo()
> R version 2.15.3 (2013-03-01)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=C LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
>
> other attached packages:
> [1] raster_2.1-25 sp_1.0-9
>
>
> loaded via a namespace (and not attached):
> [1] grid_2.15.3 lattice_0.20-15
> >
> > ##Write tile indexed MODIS file
> > #Prelims
> > nColVal<-1900
> > outFileName <- 'tmp.grd'
> >
> > #Set up raster
> > out <- raster(nrow=20000,ncol=nColVal)#modisRasterAll
> >
> > ##Get number of colums and blocksize
> > outCol <- ncol(out)
> > bs <- blockSize(out)
> >
> > ##Write to file
> > out <- writeStart(out,outFileName,overwrite=TRUE)
> > for(i in 1:bs$n){
> + out <-
> + writeValues(out,cellFromRowCol(out,bs$row[i],1) :
> cellFromRowCol(out,bs$row[i] + bs$nrow[i] -1 , ncol(out)) ,bs$row[i])
> + }
> > out <- writeStop(out)
> >
> > out2 <- raster(out)
> > out2[]<- 1:ncell(out2)
> >
> > res <- out-out2
> rgdal: version: 0.8-8, (SVN revision 463)
> Geospatial Data Abstraction Library extensions to R successfully loaded
> Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
> Path to GDAL shared files: /usr/share/gdal/1.9
> Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
> Path to PROJ.4 shared files: (autodetected)
> >
> > sum(res[])
> [1] 0
> >
> > tail(out)
> 1881 1882 1883 1884 1885 1886 1887
> 1888
> 19991 37982880 37982880 37982884 37982884 37982884 37982888 37982888
> 37982888
> 19992 37984780 37984784 37984784 37984784 37984784 37984784 37984788
> 37984788
> 19993 37986680 37986680 37986684 37986684 37986684 37986688 37986688
> 37986688
> 19994 37988580 37988584 37988584 37988584 37988584 37988584 37988588
> 37988588
> 19995 37990480 37990480 37990484 37990484 37990484 37990488 37990488
> 37990488
> 19996 37992380 37992384 37992384 37992384 37992384 37992384 37992388
> 37992388
> 19997 37994280 37994280 37994284 37994284 37994284 37994288 37994288
> 37994288
> 19998 37996180 37996184 37996184 37996184 37996184 37996184 37996188
> 37996188
> 19999 37998080 37998080 37998084 37998084 37998084 37998088 37998088
> 37998088
> 20000 37999980 37999984 37999984 37999984 37999984 37999984 37999988
> 37999988
> 1889 1890 1891 1892 1893 1894 1895
> 1896
> 19991 37982888 37982888 37982892 37982892 37982892 37982896 37982896
> 37982896
> 19992 37984788 37984792 37984792 37984792 37984792 37984792 37984796
> 37984796
> 19993 37986688 37986688 37986692 37986692 37986692 37986696 37986696
> 37986696
> 19994 37988588 37988592 37988592 37988592 37988592 37988592 37988596
> 37988596
> 19995 37990488 37990488 37990492 37990492 37990492 37990496 37990496
> 37990496
> 19996 37992388 37992392 37992392 37992392 37992392 37992392 37992396
> 37992396
> 19997 37994288 37994288 37994292 37994292 37994292 37994296 37994296
> 37994296
> 19998 37996188 37996192 37996192 37996192 37996192 37996192 37996196
> 37996196
> 19999 37998088 37998088 37998092 37998092 37998092 37998096 37998096
> 37998096
> 20000 37999988 37999992 37999992 37999992 37999992 37999992 37999996
> 37999996
> 1897 1898 1899 1900
> 19991 37982896 37982896 37982900 37982900
> 19992 37984796 37984800 37984800 37984800
> 19993 37986696 37986696 37986700 37986700
> 19994 37988596 37988600 37988600 37988600
> 19995 37990496 37990496 37990500 37990500
> 19996 37992396 37992400 37992400 37992400
> 19997 37994296 37994296 37994300 37994300
> 19998 37996196 37996200 37996200 37996200
> 19999 37998096 37998096 37998100 37998100
> 20000 37999996 38000000 38000000 38000000
> > tail(out2)
> 1881 1882 1883 1884 1885 1886 1887
> 1888
> 19991 37982881 37982882 37982883 37982884 37982885 37982886 37982887
> 37982888
> 19992 37984781 37984782 37984783 37984784 37984785 37984786 37984787
> 37984788
> 19993 37986681 37986682 37986683 37986684 37986685 37986686 37986687
> 37986688
> 19994 37988581 37988582 37988583 37988584 37988585 37988586 37988587
> 37988588
> 19995 37990481 37990482 37990483 37990484 37990485 37990486 37990487
> 37990488
> 19996 37992381 37992382 37992383 37992384 37992385 37992386 37992387
> 37992388
> 19997 37994281 37994282 37994283 37994284 37994285 37994286 37994287
> 37994288
> 19998 37996181 37996182 37996183 37996184 37996185 37996186 37996187
> 37996188
> 19999 37998081 37998082 37998083 37998084 37998085 37998086 37998087
> 37998088
> 20000 37999981 37999982 37999983 37999984 37999985 37999986 37999987
> 37999988
> 1889 1890 1891 1892 1893 1894 1895
> 1896
> 19991 37982889 37982890 37982891 37982892 37982893 37982894 37982895
> 37982896
> 19992 37984789 37984790 37984791 37984792 37984793 37984794 37984795
> 37984796
> 19993 37986689 37986690 37986691 37986692 37986693 37986694 37986695
> 37986696
> 19994 37988589 37988590 37988591 37988592 37988593 37988594 37988595
> 37988596
> 19995 37990489 37990490 37990491 37990492 37990493 37990494 37990495
> 37990496
> 19996 37992389 37992390 37992391 37992392 37992393 37992394 37992395
> 37992396
> 19997 37994289 37994290 37994291 37994292 37994293 37994294 37994295
> 37994296
> 19998 37996189 37996190 37996191 37996192 37996193 37996194 37996195
> 37996196
> 19999 37998089 37998090 37998091 37998092 37998093 37998094 37998095
> 37998096
> 20000 37999989 37999990 37999991 37999992 37999993 37999994 37999995
> 37999996
> 1897 1898 1899 1900
> 19991 37982897 37982898 37982899 37982900
> 19992 37984797 37984798 37984799 37984800
> 19993 37986697 37986698 37986699 37986700
> 19994 37988597 37988598 37988599 37988600
> 19995 37990497 37990498 37990499 37990500
> 19996 37992397 37992398 37992399 37992400
> 19997 37994297 37994298 37994299 37994300
> 19998 37996197 37996198 37996199 37996200
> 19999 37998097 37998098 37998099 37998100
> 20000 37999997 37999998 37999999 38000000
> > tail(res)
> 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894
> 19991 -1 -2 1 0 -1 2 1 0 -1 -2 1 0 -1 2
> 19992 -1 2 1 0 -1 -2 1 0 -1 2 1 0 -1 -2
> 19993 -1 -2 1 0 -1 2 1 0 -1 -2 1 0 -1 2
> 19994 -1 2 1 0 -1 -2 1 0 -1 2 1 0 -1 -2
> 19995 -1 -2 1 0 -1 2 1 0 -1 -2 1 0 -1 2
> 19996 -1 2 1 0 -1 -2 1 0 -1 2 1 0 -1 -2
> 19997 -1 -2 1 0 -1 2 1 0 -1 -2 1 0 -1 2
> 19998 -1 2 1 0 -1 -2 1 0 -1 2 1 0 -1 -2
> 19999 -1 -2 1 0 -1 2 1 0 -1 -2 1 0 -1 2
> 20000 -1 2 1 0 -1 -2 1 0 -1 2 1 0 -1 -2
> 1895 1896 1897 1898 1899 1900
> 19991 1 0 -1 -2 1 0
> 19992 1 0 -1 2 1 0
> 19993 1 0 -1 -2 1 0
> 19994 1 0 -1 2 1 0
> 19995 1 0 -1 -2 1 0
> 19996 1 0 -1 2 1 0
> 19997 1 0 -1 -2 1 0
> 19998 1 0 -1 2 1 0
> 19999 1 0 -1 -2 1 0
> 20000 1 0 -1 2 1 0
> >
> > res[8831,]
> [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [25] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [49] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [73] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [97] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [121] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [145] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [169] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [193] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
> 0 0
> [217] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [241] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [265] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [289] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [313] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [337] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [361] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [385] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [409] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
> [433] -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0 1 0 -1 0
> 1 0
>
>
> >>> Ted Rosenbaum 04/30/13 3:18 PM >>>
> Dear Matteo,
>
> Thanks for your response.
>
> You approach would work except that my computer can't handle it well since
> it loads everything into memory. That is why I used the more complicated
> approach outlined here (
> http://cran.r-project.org/web/packages/raster/vignettes/functions.pdf).
>
>
>
>
>
>
> On Tue, Apr 30, 2013 at 3:24 AM, Matteo Mattiuzzi <
> matteo.mattiuzzi at boku.ac.at> wrote:
>
> > Dear Ted,
> >
> >
> > Isn't it simply that?
> > out <- raster(nrow=20000,ncol=1900)#modisRasterAll
> > out[] <- 1:ncell(out)
> >
> >
> > tail(out)
> >
> >
> >
> >
> > >>> Ted Rosenbaum 04/29/13 11:28 PM >>>
> > Hi,
> >
> > I am trying to assign a unique ID number to each tile in a large raster
> > file.
> >
> > My code is below:
> >
> > ##Write tile indexed MODIS file
> > #Prelims
> > nColVal<-190
> > outFileName <- 'tmp.grd'
> >
> > #Set up raster
> > out <- raster(nrow=20000,ncol=nColVal)#modisRasterAll
> >
> > ##Get number of colums and blocksize
> > outCol <- ncol(out)
> > bs <- blockSize(out)
> >
> > ##Write to file
> > out <- writeStart(out,outFileName,overwrite=TRUE)
> > for(i in 1:bs$n){
> > print(i)
> > out <-
> >
> >
> writeValues(out,((outCol*(bs$row[i]-1))+1):((bs$row[i]+bs$nrows[i]-1)*outCol),bs$row[i])
> > }
> > out <- writeStop(out)
> >
> > tail(out)
> >
> > When I run this code with nColVal=190, I get the expected result from
> > tail(out): a sequence of numbers. However, when I run this code with
> > nColVal=1900, I do not and instead get
> > "37999996 38000000 38000000 38000000" for out[20000,1897:1900].
> >
> > I would like to get instead "37999997 37999998 37999999 38000000".
> >
> > Thanks for any assistance,
> >
> > Ted Rosenbaum
> > Department of Economics
> > Yale University
> >
> > [[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
> >
> >
> >
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list