# Georeferencing¶

A point on earth can be defined unambiguously in two ways:

- DATUM (lon, lat, datum)
- longitudes and latitudes are angular coordinates of a point on an ellipsoid (datum)
- PROJ (x, y, projection)
- x (eastings) and y (northings) are cartesian coordinates of a point in a map projection (the unit of x, y is usually meter)

Salem adds a third coordinate reference system (crs) to this list:

- GRID (i, j, Grid)
- on a structured grid, the (x, y) coordinates are distant of a constant (dx, dy) step. The (x, y) coordinates are therefore equivalent to a new reference frame (i, j) proportional to the projection’s (x, y) frame.

Transformations between datums and projections is handled by several tools in the python ecosystem, for example GDAL or the more lightweight pyproj, which is the tool Salem is using internally [1].

The concept of Grid added by Salem is useful when transforming data between two structured datasets, or from an unstructured dataset to a structured one.

[1] | Most datasets nowadays are defined in the WGS 84 datum, therefore the concepts of datum and projection are often interchangeable: (lon, lat) coordinates are equivalent to cartesian (x, y) coordinates in the plate carree projection. |

## Grids¶

A `Grid` is defined by a projection, a reference point in
this projection, a grid spacing and a number of grid points:

```
In [1]: import numpy as np
In [2]: import salem
In [3]: from salem import wgs84
In [4]: grid = salem.Grid(nxny=(3, 2), dxdy=(1, 1), ll_corner=(0.5, 0.5), proj=wgs84)
In [5]: x, y = grid.xy_coordinates
In [6]: x
Out[6]:
array([[ 0.5, 1.5, 2.5],
[ 0.5, 1.5, 2.5]])
In [7]: y
Out[7]:
array([[ 0.5, 0.5, 0.5],
[ 1.5, 1.5, 1.5]])
```

The default in Salem is to define the grids according to the pixels center point:

```
In [8]: smap = salem.Map(grid)
In [9]: smap.set_data(np.arange(6).reshape((2, 3)))
In [10]: lon, lat = grid.ll_coordinates
In [11]: smap.set_points(lon.flatten(), lat.flatten())
In [12]: smap.visualize(addcbar=False)
```

But with the `pixel_ref` keyword you can use another convention. For Salem,
the two conventions are identical:

```
In [13]: grid_c = salem.Grid(nxny=(3, 2), dxdy=(1, 1), ll_corner=(0, 0), proj=wgs84, pixel_ref='corner')
In [14]: assert grid_c == grid
```

While it’s good to know how grids work, most of the time grids should be inferred directly from the data files (if not, we have to implement it!):

```
In [15]: ds = salem.open_xr_dataset(salem.get_demo_file('himalaya.tif'))
In [16]: grid = ds.salem.grid
In [17]: grid.proj.srs
Out[17]: '+units=m +init=epsg:4326 '
In [18]: grid.extent
Out[18]:
[78.312513473116013,
93.937514288103202,
24.237496572347762,
31.545830286884001]
```

Grids come with several convenience functions, for example for transforming a point onto the grid coordinates:

```
In [19]: grid.transform(85, 27, crs=salem.wgs84)
Out[19]: (801.99834136842378, 544.9996059728029)
```

Or for reprojecting structured data (the xarray accessors
`transform()` method calls
`map_gridded_data()` internally):

```
In [20]: dse = salem.open_xr_dataset(salem.get_demo_file('era_interim_tibet.nc'))
In [21]: t2_era_reproj = ds.salem.transform(dse.t2m.isel(time=0))
In [22]: t2_era_reproj.salem.quick_map()
Out[22]: <salem.graphics.Map at 0x7f028bd6bbd0>
```