Map projections is a systematic transformation a location on the surface of the earth or a portion of the earth (sphere or an ellipsoid) into on a plane (flat piece of paper or computer screen).
The coordinate reference system (CRS) or spatial reference system (SRS) is a coordinate based system used to locate features on the earth.
In this section., we will cover following topics:
Data for this lesson could be found here
Geographic coordinate system (GCS)¶
Geographic coordinate systems use latitude and longitude to measure and locate feature on the globe. The GCS defines a position as a function of direction and distance from a center point of the globe, where the units of measurement are degrees. Any location on earth can be referenced by a point with longitude and latitude coordinates. For example, below figure shows a geographic coordinate system where a location is represented by the coordinate's longitude 80 degree East and latitude 55 degree North.
The equidistant lines that run east and west each have a constant latitude value called parallels. The equator is the largest circle and divides the earth in half. It is equal in distance from each of the poles, and the value of this latitude line is zero. Locations north of the equator has positive latitudes that range from 0 to +90 degrees, while locations south of the equator have negative latitudes that range from 0 to -90 degrees.
The lines that run through north and south each have a constant longitude value and form circles of the same size around the earth known as meridians. The prime meridian is the line of longitude that defines the origin (zero degrees) for longitude coordinates. One of the most commonly used prime meridian locations is the line that passes through Greenwich, England. Locations east of the prime meridian up to its antipodal meridian (the continuation of the prime meridian on the other side of the globe) have positive longitudes ranging from 0 to +180 degrees. Locations west of the prime meridian has negative longitudes ranging from 0 to -180 degrees.
The latitude and longitude lines cover the globe to form a grid that know as a graticule. The point of origin of the graticule is (0,0), where the equator and the prime meridian intersect. The equator is the only place on the graticule where the linear distance corresponding to one degree latitude is approximately equal the distance corresponding to one degree longitude. Because the longitude lines converge at the poles, the distance between two meridians is different at every parallel.
The World Geodetic system 84¶
The most recent geographic coordinate system is the World Geodetic system 84 also known as WGS 1984 or EPSG:4326 (EPSG- European Petroleum Survey Group). It consists of a standard coordinate system, spheroidal reference (the datum or reference ellipsoid) and raw altitude.
Projected coordinate system¶
A projected coordinate system provide mechanisms to project maps of the earth's spherical surface onto a two-dimensional Cartesian coordinate (x,y coordinates) plane. Projected coordinate systems are referred to as map projections. This approach is useful where accurate distance, angle, and area measurements are needed. The term 'projection' is often used interchangeably with projected coordinate systems.
Commonly use projected coordinate systems include:
Universal Transverse Mercator
Albers Equal Area Conic
Universal Transverse Mercator¶
The most widely used two-dimensional Cartesian coordinate system is the Universal Transverse Mercator (UTM) system which represents a horizontal position on the globe and can be used to identify positions without having to know their vertical location on the 'y' axis. The UTM system is not a single map projection. It represents the earth as sixty different zones, each composed of six-degree longitudinal bands, with a secant transverse Mercator projection in each. The UTM system divides the Earth between 80 S and 84 N latitude into 60 zones, each 6 of longitude in width. Zone 1 covers longitude 180 to 174 W; zone numbering increases eastward to zone 60, which covers longitude 174 to 180.
The Universal Transverse Mercator grid that covers the conterminous 48 United States comprises 10 zones, from Zone 10 on the west coast through Zone 19 in New England USGS, 2001.The New York State falls in zones 17N (EPSG 26917) and zone 18N (EPSG: 26918) with North American Datum 1983 (NAD83).
Albers Equal Area Conic¶
Albers Equal Area Conic projection system uses two standard parallels to reduce the distortion of shape in the region between the standard parallels. This projection is best suited for area extending in an east-to-west orientation like the conterminous United States. Used for the conterminous United States, normally using 29-30 and 45-30 as the two standard parallels. For this projection, the maximum scale distortion for the 48 states of the conterminous United States is 1.25 percent.
Load Libraries¶
import os
from glob import glob
import geopandas as gpd
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import geopandas as gpd
import rasterio as rio
from rasterio.plot import show
import rasterio.warp
import rasterio.shutil
from rasterio.warp import calculate_default_transform, reproject, Resampling
Set working directory¶
path= "C:/Users/zahmed2/Documents/spatial_data_processing_and_analysis_python/Lesson_03_map_pojection_coordinate_reference_system/"
os.chdir(path)
print("Current Working Directory " , os.getcwd())
Coordinate Reference System in Python¶
Vector Data¶
# Import NY Ctate County Shape file
fp_poly_gcs="Data_03//NY_County_GCS.shp"
poly_gcs= gpd.read_file(fp_poly_gcs)
Check CRS¶
View the CRS of each layer using .crs attribute of geopandas dataframes (e.g. dataframename.crs).
print(poly_gcs.crs)
We noticed that CRS of NYS county shape has not been be defined yet. So you need to define its current CRS (WGS 1984 or EPSG:4326) before you do any further analyses like re-projection etc. Yoy would set the projection by assigning the WGS84 latitude-longitude CRS (EPSG:4326) to the crs attribute:
poly_gcs.crs = {'init' :'epsg:4326'}
print(poly_gcs.crs)
ax1=poly_gcs.plot()
# add a title
ax1.set_title('WGS84 (lat/lon)')
Re-pojection¶
We will transform county shape file from WGS 1984 to Albers Equal Area Conic NAD1983 (EPSG:5070). You could reproject the data from one CRS to another CRS using the geopandas method (i.e. function belonging to geopandas objects): to_crs(). Function .to_crs() will only work if your original spatial object has a CRS assigned to it AND if that CRS is the correct CRS for that data!
GeoPandas uses pyproj for the projections (requirements) and pyproj does not know the crs EPSG:5070 (does not exist) nor the ESRI:10203. But you can try the Proj4 string of the ESRI:10203 projection
Polygoan¶
# Define EPSG:5070 projection
new_crs={'init': 'epsg:5070'}
poly_proj=poly_gcs.to_crs(new_crs)
ax2=poly_proj.plot()
# add a title
ax2.set_title('WGS84 (lat/lon)')
# Make subplots that are next to each other
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 8))
# Plot the data in WGS84 CRS
poly_gcs.plot(ax=ax1, facecolor='gray');
# Add title
ax1.set_title("WGS84");
# Plot the one with ETRS-LAEA projection
poly_proj.plot(ax=ax2, facecolor='blue');
# Add title
ax2.set_title("Albers Equal Area Conic");
# Remove empty white space around the plot
plt.tight_layout()
Point shape file¶
# Load point shape file
fp_point_gcs="Data_03/GP_Point_GCS.shp"
point_gcs= gpd.read_file(fp_point_gcs)
# Define projection
point_gcs.crs = {'init' :'epsg:4326'}
print(point_gcs.crs)
# plot map
ax3=point_gcs.plot()
ax3.set_title('WGS84 (lat/lon)')
# Change to EPSG:5070 projection
point_proj=point_gcs.to_crs(new_crs)
ax4=point_proj.plot()
ax4.set_title('Albers Equal Area Conic)')
## Shave projected file
point_proj.to_file("GP_Point_PROJ.shp")
Raster Data¶
This exercise you will learn how to projection transformation will be done of raster data. In this exercise, we will use SRTM 90 digital Elevation Model New York State which was downloaded from CGIAR-CSI. We will re-project it from WSG84 coordinate to Albers Equal Area Conic NAD83 projection system.
Projection of a Single Raster¶
dem_path = "Data_03/Onondaga_DEM.tif"
dem=rio.open(dem_path)
show(dem)
Read raster file properties¶
# Check projection
dem.crs
# Check projection with GDAL
import gdal
!gdalinfo "Data_03/Onondaga_DEM.tif"
# All Metadata for the whole raster dataset
dem.meta
First, we will re-project one tile of DEM raster data from WGS 1984 CRS to Albers Equal Area Conic projection system. In Python, re-projection of raster little complicated. The code below reprojects between two arrays, using no pre-existing GIS datasets. The rasterio.warp.reproject() has two positional arguments: source and destination. The remaining keyword arguments parameterize the reprojection transform.
Reprojecting a GeoTIFF dataset from one coordinate reference system is a common use case. Rasterio provides a few utilities to make this even easier:
transform_bounds() transforms the bounding coordinates of the source raster to the target coordinate reference system, densifiying points along the edges to account for non-linear transformations of the edges.
calculate_default_transform() transforms bounds to target coordinate system, calculates resolution if not provided, and returns destination transform and dimensions.
inputTIF = 'Data_03/Onondaga_DEM.tif'
projTIF = 'Data_03/Onondaga_DEM_PROJ.tif'
dst_crs = new_crs
with rio.open(inputTIF) as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rio.open(projTIF, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rio.band(src, i),
destination=rio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
with rio.open(projTIF) as src:
rasterio.plot.show(src, title='Albers Equal Area Conic', cmap='RdYlGn')
Below reproject_et function wrapps the rasterio reproject code with some addtional arguments that allow to define input path, the output path, and the new CRS components. These three argument only need to change we are going to cahnge project other raster layers. (Source)
def reproject_et(inpath, outpath, new_crs):
dst_crs = new_crs # CRS for web meractor
with rio.open(inpath) as src:
transform, width, height = calculate_default_transform(
src.crs, dst_crs, src.width, src.height, *src.bounds)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rio.open(outpath, 'w', **kwargs) as dst:
for i in range(1, src.count + 1):
reproject(
source=rio.band(src, i),
destination=rio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=transform,
dst_crs=dst_crs,
resampling=Resampling.nearest)
reproject_et(inpath = 'Data_03/Onondaga_DEM.tif',
outpath = 'Data_03/Onondaga_DEM_PROJ_02.tif',
new_crs = "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs")
# Check projection with GDAL
import gdal
!gdalinfo "Data_03/Onondaga_DEM_PROJ_02.tif"