From Degrees To Meters: Transforming Geographic Data To Projected Crs With Qgis

Understanding Coordinate Reference Systems

A coordinate reference system (CRS) defines how a two- or three-dimensional space is divided into coordinates that allow the positions of points in that space to be uniquely defined. Geographic coordinate systems use latitude and longitude in angular units like degrees while projected coordinate systems use linear units like meters.

The most common global geographic CRS is the World Geodetic System 1984 (WGS84) which defines coordinates on an ellipsoidal model of Earth. However, analysis operations like distance and area calculations are more accurately performed on projected planar surfaces. Transforming data from geographic to projected CRS enables such processing while retaining the ability to integrate with data referenced to global spherical coordinates.

Preparing Data in QGIS

QGIS supports hundreds of standard coordinate systems and seamless on-the-fly transformation between them. Loading spatial data, identifying native coordinate systems, and setting an appropriate project CRS are key initial steps.

Vector data such as ESRI shapefiles and GeoPackage tables with geometry columns can be loaded via Layer > Add Layer > Add Vector Layer. File database sources like PostGIS are also available. Review a layer’s properties to verify its native CRS in the Information section. If undefined, the project CRS will be assumed.

Raster data from formats like GeoTIFF and JPEG2000 can be added with Layer > Add Layer > Add Raster Layer. Check raster properties and note whether pixel values represent locations implicitly associated with a particular CRS or if georeferencing information is defined in an included world file.

For global data integration, set Project Properties > CRS to a standard geographic system like WGS 84 (EPSG:4326). Enable on-the-fly reprojection for layers with differing native CRS. QGIS will transform coordinates to the set project system without altering original data.

Transforming Vector Data

To explicitly convert vector data to another CRS, use the Layer > Save As… menu. This creates a physical copy of the data with coordinates changed to the selected target system. Common projected targets are UTM and USA Contiguous Albers Equal Area Conic to optimize metric distance and area analysis tasks.

Alternatively, reprojection on-the-fly allows temporary analysis views of data in varying CRS without duplication. Enable this mode in Project Properties under the CRS section. Make sure “Enable ‘on the fly’ CRS transformation” is checked. The project CRS will be used for display and measurement.

For example, setting a UTM projection and overlaying originally lat/long wildlife tracking data facilitates easier distance comparison. The underlying data remains unaltered, reverting when disabling on-the-fly transformation.

Transforming Raster Data

To explicitly reproject raster data, use the Raster > Projections > Warp (Reproject) tool. This resamples cell values from the source CRS into a target CRS based on user settings.

Select an appropriate resampling algorithm based on your data and analysis needs.Nearest neighbour conversion is faster but can lead to distortion while cubic interpolation is smoother but slower. For continuous data like elevation models, bilinear or cubic is recommended. Categorical or integer-based data like land cover grids can use nearest neighbour.

Additionally, clipping rasters to a rectangular extent that encloses the area of interest helps avoid excess null data that increases file sizes without additional information value.

Checking Transformed Data Accuracy

After reprojecting geographic data, inspect visually against basemaps or orthophoto imagery to ensure proper alignment. Source control points like surveyed locations can also be checked. For example, overlay a reprojected vector stream layer on high resolution aerial photos. Do features align with true locations or is there noticeable offset?

For point layers, inspect locations of a systematic subset of coordinates in the attribute table before and after transformation to verify latitude/longitude or easting/northing values convert appropriately without introducing unexpected shifts or distortion.

If alignment issues occur, review projection metadata and parameters to identify and correct discrepancies. Ensure accurate selection of appropriate projection datums and zones to suit the geographic area covered by your data.

Example Code for Vector Transformation

Here is Python code to batch reproject a folder of ESRI shapefiles from WGS 84 to UTM Zone 15N based on the WGS84 datum:

import os
import fiona
from fiona.crs import from_epsg
from pyproj import transform, Proj

input_crs = from_epsg(4326)
output_crs = from_epsg(32615)

transform_function = lambda x, y: transform(Proj(init='epsg:4326'), Proj(init='epsg:32615'), x, y) 

for filename in os.listdir('path/to/shapefiles'):
    name, ext = os.path.splitext(filename)
    if ext.lower() == '.shp':

        with fiona.open('path/to/shapefiles/' + filename) as source:

            output_path = 'output/reprojected/' + name + '_32615.shp'        
            with fiona.open(output_path, 'w', driver='ESRI Shapefile', crs=output_crs, schema=source.schema) as output:
            
                for feature in source:
                
                    new_coords = [transform_function(lon, lat) for lon, lat in feature['geometry']['coordinates']]
                    feature['geometry']['coordinates'] = new_coords

                    output.write(feature)

print('Transformation complete')

Example Code for Raster Transformation

Here is Python code using GDAL to batch reproject and clip a folder of GeoTIFF rasters from WGS 84 to UTM Zone 10N based on the NAD83 datum:

import os
from osgeo import gdal, ogr, osr

input_crs = osr.SpatialReference()
input_crs.ImportFromEPSG(4326)

output_crs = osr.SpatialReference() 
output_crs.ImportFromEPSG(26910)

envelope = ogr.CreateGeometryFromWkt('POLYGON ((500000 4500000, 500000 4510000, 600000 4510000, 600000 4500000, 500000 4500000))')

for filename in os.listdir('/rasters'):
    name, ext = os.path.splitext(filename)
    if ext == '.tif':

        src = gdal.Open('/rasters/' + filename)
       
        output = '/output/' + name + '_26910.tif'
       
        warp_options = gdal.WarpOptions(dstSRS = output_crs, 
                                        outputBounds = envelope.GetEnvelope(),  
                                        xRes = 30, yRes = 30)
                                        
        result = gdal.Warp(output, src, options = warp_options)
        
        src = None
        result = None
        
print('Processing complete')

Leave a Reply

Your email address will not be published. Required fields are marked *