# Title 
# -- Read SIMPLE-G-US-2010-v2 data from a text file

# Prepared by 
# -- Iman Haqiqi, ihaqiqi@purdue.edu
# -- Global to Local Analysis of Systems Sustainability (GLASS)
# -- Global Trade Analysis Project (GTAP)
# -- Department of Agricultural Economics, Purdue University

# load required libraries
library(terra)

# read the text file which includes longitude and latitude
df = read.csv("SIMPLEG_QCROP_QLAND_5min_US_2010_v2.csv", sep=",", header=T)

# reproduce the coordinates for 5-arc min to avoid decimal precision issues 
df$x = round(df$x * 120) / 120
df$y = round(df$y * 120) / 120

head(df)
#           x        y  FIPS SubRegion QCROPgl_irr QCROPgl_rfd QLANDgl_irr QLANDgl_rfd
# 1 -124.4583 42.79167 41015         7       0.305       0.101       0.109        0.01
# 2 -124.4583 42.87500 41015         7       0.794       0.101       0.284        0.01
# 3 -124.3750 42.45833 41015         7       0.202       0.101       0.010        0.01
# 4 -124.3750 42.95833 41011         7       0.202       0.101       0.010        0.01
# 5 -124.3750 43.04167 41011         7       0.680       0.101       0.281        0.01
# 6 -124.3750 43.12500 41011         7       0.878       0.101       0.363        0.01


# convert the dataframe to a raster stack
r = rast(df, type="xyz")

# save as raster files file
writeRaster(r$QCROPgl_irr, "QCROPgl_irr_kMT_5min_2010_v2.tif")
writeRaster(r$QCROPgl_rfd, "QCROPgl_rfd_kMT_5min_2010_v2.tif")
writeRaster(r$QLANDgl_irr, "QLANDgl_irr_kHA_5min_2010_v2.tif")
writeRaster(r$QLANDgl_rfd, "QLANDgl_rfd_kHA_5min_2010_v2.tif")

writeCDF(r$QCROPgl_irr, 
         "QCROPgl_irr_kMT_5min_2010_v2.nc", 
         varname="QCROPgl_irr", unit="1000 Metric Ton", 
         longname="corn-equivalent crop production circa 2010, irrigated")

writeCDF(r$QCROPgl_rfd, 
         "QCROPgl_rfd_kMT_5min_2010_v2.nc", 
         varname="QCROPgl_rfd", unit="1000 Metric Ton", 
         longname="corn-equivalent crop production circa 2010, rainfed")

writeCDF(r$QLANDgl_irr, 
         "QLANDgl_irr_kHA_5min_2010_v2.nc", 
         varname="QLANDgl_irr", unit="1000 Hectare", 
         longname="cropland area circa 2010, irrigated")

writeCDF(r$QLANDgl_rfd, 
         "QLANDgl_rfd_kHA_5min_2010_v2.nc", 
         varname="QLANDgl_rfd", unit="1000 Hectare", 
         longname="cropland area circa 2010, rainfed")