library(MODISTools)
library(raster)
library(gdalUtils)
library(RCurl)
setwd("/home/slavko/landcover/turska")
wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

#;MOD13A3 <- "https://e4ftl01.cr.usgs.gov/MODV6_Cmp_B/MOLT/MOD13A3.006/2018.01.01/"
#;items <- strsplit(getURL(MOD13A3), "\n")[[1]]


sds12 <- get_subdatasets("MOD13A3.A2018001.h19v06.006.2018033223418.hdf")
sds11 <- get_subdatasets("MOD13A3.A2018001.h19v05.006.2018033223835.hdf")
sds10 <- get_subdatasets("MOD13A3.A2018001.h19v04.006.2018033223911.hdf")
sds9 <- get_subdatasets("MOD13A3.A2018001.h22v06.006.2018033223754.hdf")
sds8 <- get_subdatasets("MOD13A3.A2018001.h22v05.006.2018033224426.hdf")
sds7 <- get_subdatasets("MOD13A3.A2018001.h22v04.006.2018033224431.hdf")
sds6 <- get_subdatasets("MOD13A3.A2018001.h21v06.006.2018033223840.hdf")
sds5 <- get_subdatasets("MOD13A3.A2018001.h21v05.006.2018033223926.hdf")
sds4 <- get_subdatasets("MOD13A3.A2018001.h21v04.006.2018033224548.hdf")
sds3 <- get_subdatasets("MOD13A3.A2018001.h20v06.006.2018033224027.hdf")
sds2 <- get_subdatasets("MOD13A3.A2018001.h20v05.006.2018033224657.hdf")
sds1 <- get_subdatasets("MOD13A3.A2018001.h20v04.006.2018033224629.hdf")

gdal_translate(sds12[1],dst_dataset = "MOD13A3.A2018001.h19v06.tif")
gdal_translate(sds11[1],dst_dataset = "MOD13A3.A2018001.h19v05.tif")
gdal_translate(sds10[1],dst_dataset = "MOD13A3.A2018001.h19v04.tif")
gdal_translate(sds9[1],dst_dataset = "MOD13A3.A2018001.h22v06.tif")
gdal_translate(sds8[1],dst_dataset = "MOD13A3.A2018001.h22v05.tif")
gdal_translate(sds7[1],dst_dataset = "MOD13A3.A2018001.h22v04.tif")
gdal_translate(sds6[1],dst_dataset = "MOD13A3.A2018001.h21v06.tif")
gdal_translate(sds5[1],dst_dataset = "MOD13A3.A2018001.h21v05.tif")
gdal_translate(sds4[1],dst_dataset = "MOD13A3.A2018001.h21v04.tif")
gdal_translate(sds3[1],dst_dataset = "MOD13A3.A2018001.h20v06.tif")
gdal_translate(sds2[1],dst_dataset = "MOD13A3.A2018001.h20v05.tif")
gdal_translate(sds1[1],dst_dataset = "MOD13A3.A2018001.h20v04.tif")

gdalbuildvrt (output.vrt="MOD13A3h.vrt",gdalfile=c("MOD13A3.A2018001.h22v06.tif",
                                                   "MOD13A3.A2018001.h22v05.tif",
                                                   "MOD13A3.A2018001.h22v04.tif",
                                                   "MOD13A3.A2018001.h21v06.tif",
                                                   "MOD13A3.A2018001.h21v05.tif",
                                                   "MOD13A3.A2018001.h21v04.tif",
                                                   "MOD13A3.A2018001.h20v06.tif",
                                                   "MOD13A3.A2018001.h20v05.tif",
                                                   "MOD13A3.A2018001.h20v04.tif",
                                                   "MOD13A3.A2018001.h19v04.tif",
                                                   "MOD13A3.A2018001.h19v05.tif",
                                                   "MOD13A3.A2018001.h19v06.tif"))

gdalinfo("MOD13A3h.vrt")


#mask1a <- crop(mask, extent(x1,x2,y1,y2))
#soil1a <- crop(soil, extent(x1,x2,y1,y2))

gdalwarp("MOD13A3h.vrt",dstfile="MOD13A3wgs.tif",t_srs=wgs1984,output_Raster=TRUE,overwrite=TRUE,verbose=TRUE,r="near")
gdal_translate("MOD13A3wgs.tif","MOD13A3wgs.xyz",of="XYZ",overwrite=TRUE)
gdalinfo("MOD13A3wgs.tif")
#dat <- read.csv("MOD13A3wgs.xyz", header = FALSE)


r <- raster("MOD13A3wgs.tif")
#r[] <- runif(ncell(r))
#coordinates
#coords <- xyFromCell(r,1:ncell(r))

#create list
#xyzlist <- list(x=coords[,'x'],y=coords[,'y'],z=as.matrix(r))

zz=t(as.matrix(r))

xlen=as.integer(ncol(zz))
ylen=as.integer(nrow(zz))
aa=c(xlen,ylen)
tii=as.numeric(zz)
write(aa,file='iwc12.txt',ncolumns=2,sep=',')
write(tii,file='iwc12.txt',ncolumns=10,sep=',',append = TRUE)
stop
#fw <- file("testbin", "wb")

writeBin(con=fw,object=as.integer(nrow(zz)))
writeBin(con=fw,object=as.vector(as.numeric(zz)))

xyz <- rasterToPoints(r)



ax=list(x=coords[,'x'])
x=t(as.matrix(ax))
y=t(as.matrix(r))
z=t(as.matrix(r))


gdal_translate("MOD13A3h29v12wgs.tif","MOD13A3h29v12wgs.xyz",of="XYZ",overwrite=TRUE)
gdal_translate("MOD13A3h30v12wgs.tif","MOD13A3h30v12wgs.xyz",of="XYZ",overwrite=TRUE)
gdal_translate("MOD13A3h30v12.tif","MOD13A3h30v12.xyz",of="XYZ",overwrite=TRUE)
# Load and plot the new .tif

rast <- raster("/home/slavko/landcover/australija/NPP2000.tif")
gdal_translate("/home/slavko/landcover/australija/NPP2000.tif","/home/slavko/landcover/australija/NPP2000.xyz",of="XYZ")

plot(rast)
wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
gdalwarp("NPP2000.tif",dstfile="NPP2000wgs.tif",t_srs=wgs1984,output_Raster=TRUE,overwrite=TRUE,verbose=TRUE,r="near")


projection(rast) <- wgs1984
#crs(rast) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" 
plot(rast)