---
title: "Transforming coordinates from crs to crs"
output: md_document
---

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(
  echo = TRUE,
  fig.path = "images/",
  warning = FALSE, message = FALSE
)
```

# Transforming satellite data from one map projection to another
> history \| Create July 2023 \| Updated August 2023<br/>

## Background

Map projections try to portray the earth's spherical surface on a flat surface. A coordinate reference system (CRS) defines how the two-dimensional, projected map relates to real places on the earth. Which map projection and CRS to use depends on the region in which you are working and the analysis you will do.

<a href="https://polarwatch.noaa.gov/" target="_blank">NOAA PolarWatch</a> distributes gridded and tabular oceanographic data for polar regions.  

Most of the satellite data in the PolarWatch data catalog use a projection based on a Geographic Coordinate Reference System, where the X and Y coordinates are longitude and latitude, respectively. Geographical coordinates work well in many parts of the globe, but within polar regions, features tend to be very distorted. A Polar Stereographic projection often is a better choice for polar regions. For example, the NSIDC's Polar Stereographic Projections, which were developed to optimize mapping of sea ice, have only a 6% distortion of the grid at the poles and no distortion at 70º, a latitude close to where the marginal ice zones occur. The <a href="https://polarwatch.noaa.gov/catalog/ice-sq-sh-nsidc-cdr-v4/preview/?dataset=daily&var=cdr_seaice_conc">NOAA NSIDC Sea Ice Concentration Climate Data Record</a> dataset, for example, is in a polar stereographic projection. 

When working with satellite datasets with a mix of map projections, it is often necessary to transform all the data to a common projection. 
In this exercise, we will learn to transform coordinates from one projection to another.

## Objective
In this tutorial, we will learn to transform dataset coordinates from one projection to another.

## This tutorial demonstrates the following techniques
- Downloading and saving a netCDF file from PolarWatch ERDDAP data server
- Accessing satellite data and metadata in polar stereographic projection 
- Convert the netcdf data into a dataframe
- Transforming coordinates using EPSG codes
- Mapping data using the transformed coordinates

## Dataset used
    
*Sea Ice Concentration, NOAA/NSIDC Climate Data Record V4, Southern Hemisphere, 25km, Science Quality, 1978-2022, Monthly*
The sea ice concentration (SIC) dataset used in this exercise is produced by NOAA NSIDC from passive microwave sensors as part of the Climate Data Record. It is a science quality dataset of monthly averages that extends from 1978-2022. SIC is reported as the fraction (0 to 1) of each grid cell that is covered by ice. The data are mapped in the Southern Polar Stereographic projection (EPSG:3031). The resolution is 25km, meaning each pixel in this data set represents a value that covers a 25km by 25km area. The dataset can be downloaded directly from the PolarWatch ERDDAP at the following link:  https://polarwatch.noaa.gov/erddap/griddap/nsidcG02202v4shmday


### Install required packages

This code block will check if required packages are installed, and will install missing packages.

```{r setup_packages, message=FALSE, warning=FALSE}

# Function to check if pkgs are installed, install missing pkgs, and load
pkgTest <- function(x)
{
  if (!require(x,character.only = TRUE))
  {
    install.packages(x,dep=TRUE,repos='http://cran.us.r-project.org')
    if(!require(x,character.only = TRUE)) stop(x, " :Package not found")
  }
}

list.of.packages <- c("ncdf4" , "sf", "ggplot2","scales", "RColorBrewer")

# create list of installed packages
pkges = installed.packages()[,"Package"]
for (pk in list.of.packages) {
  pkgTest(pk)
}

```

### Get data from ERDDAP 
```{r}
# download the sea ice data NetCDF file
url <- "https://polarwatch.noaa.gov/erddap/griddap/nsidcG02202v4shmday.nc?cdr_seaice_conc_monthly[(2022-12-01T00:00:00Z):1:(2022-12-01T00:00:00Z)][(4350000.0):1:(-3950000.0)][(-3950000.0):1:(3950000.0)]"

sic <- download.file(url, destfile="sic.nc", mode='wb')

# file open
ds <- nc_open('sic.nc')

# print metadata
print(ds)
```


```{r}
# get data into r variables 
xgrid <- ncvar_get(ds, "xgrid")
ygrid <- ncvar_get(ds, "ygrid")
sic <- ncvar_get(ds, "cdr_seaice_conc_monthly")  #lat and lon
dim(sic)


# close 
nc_close(ds)
```
### Map the polar stereographic projected data 

```{r}

# create dataframe
sicd <- expand.grid(xgrid=xgrid, ygrid=ygrid)
sicd$sic <- array(sic, dim(xgrid)*dim(ygrid))

# exclude fillvalue
sicd$sic[sicd$sic > 2] <- NA 

# map sea ice concentration
ggplot(data = sicd, aes(x = xgrid, y = ygrid, fill=sic) ) + 
       geom_tile() + 
       coord_fixed(ratio = 1) + 
       scale_y_continuous(labels=comma) + 
       scale_x_continuous(labels=comma) +
       scale_fill_gradientn(colours=rev(brewer.pal(n = 3, name = "Blues")),na.value="tan")+ggtitle("Sea Ice Concentration in Polar Steregraphic projection")
```


## Transforming from CRS to CRS

When transforming from one CRS to another, it is important to inspect CRS definitions and the transformation function for proper transformation.
We will transform from CRS EPSG: 3031 (NSIDC Polar Stereographic South) to EPSG: 4326 (geographic coordinate system) 

```{r}

xy_sf <- st_as_sf(sicd, coords=c("xgrid", "ygrid"), crs=3031)

# st_transform output order: lat , lon
latlon_sf <- st_transform(xy_sf, crs=4326)

latlon_df <- as.data.frame(st_coordinates(latlon_sf))
names(latlon_df) <- c("Lat", "Lon")
seaiceconc <- sicd$sic
# create dataframe and add names
sicdf_latlon <- data.frame(cbind(latlon_df,seaiceconc))
head(na.omit(sicdf_latlon), 5)


```
### Plot data with new coordinates on a global map

```{r}
sic_sf <- st_as_sf(sicdf_latlon, coords = c('Lat', 'Lon'))
st_crs(sic_sf) = "4326"

# plot sea ice concentration on an unprojected map
ggplot(sic_sf) + geom_sf(aes(color = seaiceconc)) + labs(title = "Sea Ice Concentration of Antarctica")



```

Because the data were in the polar stereographic projection, mapping the data onto a different projection (global map projection) above, doesn't make the data fit well on the map.

## References

- <a href="https://polarwatch.noaa.gov/catalog/ice-sq-nh-nsidc-cdr-v4/preview/?dataset=daily&var=cdr_seaice_conc&time_min=2022-05-31T00:00:00Z&time_max=2022-05-31T00:00:00Z&proj=epsg3413&colorBar=KT_ice|||0|1|" target="_blank">NOAA PolarWatch Data Product Page (download, preview)</a>
- <a href="https://nsidc.org/data/g02202/versions/4" target="_blank">NSIDC Data Product Description</a>
- <a href="https://nsidc.org/sites/default/files/g02202-v004-userguide_1_1.pdf" target="_blank">NSIDC Data Product User Guide (pdf)</a>