I am struggling to automate this.
### script to retrieve first 3 days of data from Oct 2011 ###
require(rgdal)
yyyy<- "2011"
mm<- "10"
for(dd in 1:3){
file <- paste("nt_",yyyy,mm,sprintf("%2.2d", dd),"_f17_nrt_n.bin",
sep="")
data.bin <- paste("
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/
", file, sep="")
#Open binary connection to the file
bincon<- file(data.bin, 'rb')
#Here I tried to pass bincon to the .vrt
ice<-readGDAL("nsidcBin.vrt")
image(ice)
close(bincon)
}
This fails, because I am not able to pass the binary reading connection
(bincon) to the vrt file
I tried to simply write the .vrt file to have the binary data explicitly
under in the SourceFilename argument
<SourceFilename relativetoVRT="0">
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/nt_20111011_f17_nrt_n.bin
</SourceFilename>
but this fails with the following error.
in .local(.Object, ...) :
Unable to open
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/nt_20111011_f17_nrt_n.bin
.
No error
Would I be better off to read the binary file in as a matrix and coerce it
into a spatial pixel dataframe from there?
require(sp)
# make loop here
yyyy<- "2011"
mm<- "10"
for(dd in 1:3){
file <- paste("nt_",yyyy,mm,sprintf("%2.2d", dd),"_f17_nrt_n.bin",
sep="")
data.bin <- paste("
ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time/north/
", file, sep="")
bincon<- file(data.bin, 'rb')
result <- matrix(0,448,304) #built result matrix with 448 rows
and 304 columns
# step past the first 300 bytes
for(i in 1:300){readBin(bincon,"raw",1)}
# step through the 448 rows of 304 columns
for(i in 1:448){
for(j in 1:304){
x <- as.integer(readBin(bincon,"raw",1))
if(x<251) result[i,j]<- x/250
}
}
close(bincon)
image(result)
}
>From here I would need to build the spatial pixel data frame from the
matrix...
Tony Fischbach, Wildlife Biologist
Walrus Research Program
Alaska Science Center
U.S. Geological Survey
4210 University Drive
Anchorage, AK 99508-4650
Phone: 907 786-7145
FAX: 907 786-7150
[email protected]
http://alaska.usgs.gov/science/biology/walrus
From:
Barry Rowlingson <[email protected]>
To:
Michael Sumner <[email protected]>
Cc:
Anthony Fischbach <[email protected]>, [email protected]
Date:
10/12/2011 10:54 AM
Subject:
Re: [R-sig-Geo] Reading National Snow and Ice Data Center binary files
Sent by:
[email protected]
On Wed, Oct 12, 2011 at 11:08 AM, Barry Rowlingson
<[email protected]> wrote:
> Annoyingly I don't think gdal's raw data vrt system can find out the
> number of pixels across/down the raster - you have to code them into
> the .vrt file. I reckon with a bit of unix (or in extremis, C code)
> craftiness it would be possible to get the info from the header of the
> .bin file and write the relevant .vrt file...
The data format details tells is the files only come in a limited
number of sizes. I grabbed a .bin file and wrote this .vrt file to sit
next to it:
<VRTDataset rasterXSize="304" rasterYSize="448">
<VRTRasterBand dataType="Byte" band="1" subClass="VRTRawRasterBand">
<SourceFilename
relativetoVRT="1">nt_20111010_f17_nrt_n.bin</SourceFilename>
<ImageOffset>300</ImageOffset>
<PixelOffset>1</PixelOffset>
<LineOffset>304</LineOffset>
</VRTRasterBand>
</VRTDataset>
Then in R:
require(rgdal)
ice=readGDAL("nt_20111010_f17_nrt_n.vrt")
image(ice)
and there's Greenland, and various other bits of rapidly melting ice...
Hope this helps - the important bits of the .vrt file are the 300
offset the raster X and Y size and LineOffset(which =XSize) and
obviously the filename...
I'd write a script to create these .vrt files and then convert
everything to a standard format like GeoTiff... Note you'll have to
add the projection info on manually... etc etc etc.
Barry
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
