library(ncdf4)
library(raster)
library(sf)
library(dplyr)


# Define the bounding box for the US (approximate)
minlon = -130
maxlon = -64
minlat = 24
maxlat = 52
bbox <- c(xmin = minlon, xmax = maxlon, ymin = minlat, ymax = maxlat)
us_extent <- extent(minlon, maxlon, minlat, maxlat)

# Read the netCDF file
nc_file <- "/Users/zac/Data/PM25_Analysis/WashU/V6GL02.02.CNNPM25.Global.202201-202212.nc"
nc <- nc_open(nc_file)

# Extract latitude, longitude, and concentration data
lat <- ncvar_get(nc, "lat")
lon <- ncvar_get(nc, "lon")
pm25 <- ncvar_get(nc, "PM25")

# Close the netCDF file
nc_close(nc)

# Set negative PM2.5 values to NA
pm25[pm25 < 0] <- NA

# Flip the PM2.5 data along the latitude dimension
pm25 <- pm25[, ncol(pm25):1]

# Create a raster object for PM2.5 data
pm25_raster <- raster(t(pm25), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat))
projection(pm25_raster) <- CRS("+proj=longlat +datum=WGS84")

# Crop the raster to the US extent
r <- crop(pm25_raster, us_extent)

#nc2raster = stack(input_nc,varname = varname)

output = '/Users/zac/Data/PM25_Analysis/WashU/6GL02.02.CNNPM25.2022.tif'
writeRaster(r,output,format = 'GTiff',overwrite = TRUE)
#Change the output path to save Geotif data.