⚠️ Our 0.1 release refactored several early-development functions for long-term stability, to update your code see here. ⚠️
Future changes will come with deprecation warnings! 🙂

Reproject a vector

Reproject a vector#

This example demonstrates the reprojection of a vector using geoutils.Vector.reproject().

We open a raster and vector.

import geoutils as gu

filename_rast = gu.examples.get_path("everest_landsat_b4_cropped")
filename_vect = gu.examples.get_path("everest_rgi_outlines")
rast = gu.Raster(filename_rast)
vect = gu.Vector(filename_vect)

The two objects are in different projections.

Driver:               GTiff
Opened from file:     /home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/stable/examples/data/Everest_Landsat/LE71400412000304SGS00_B4_cropped.tif
Filename:             /home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/stable/examples/data/Everest_Landsat/LE71400412000304SGS00_B4_cropped.tif
Loaded?               False
Modified since load?  False
Grid size:                 492, 315
Number of bands:      1
Data types:           uint8
Coordinate system:    ['EPSG:32645']
Nodata value:         None
Pixel interpretation: Area
Pixel size:           30.0, 30.0
Upper left corner:    483430.0, 3093260.0
Lower right corner:   498190.0, 3102710.0


"Filename:           /home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/stable/examples/data/Everest_Landsat/15_rgi60_glacier_outlines.gpkg \nCoordinate System:  EPSG:4326\nExtent:             [86.70393205700003, 27.847778695000045, 87.11409458600008, 28.14549759700003] \nNumber of features: 86 \nAttributes:         ['RGIId', 'GLIMSId', 'BgnDate', 'EndDate', 'CenLon', 'CenLat', 'O1Region', 'O2Region', 'Area', 'Zmin', 'Zmax', 'Zmed', 'Slope', 'Aspect', 'Lmax', 'Status', 'Connect', 'Form', 'TermType', 'Surging', 'Linkages', 'Name', 'geometry']"

Let’s plot the two in their original projection.

rast.plot(cmap="Greys_r")
vect.plot(ax="new", fc="none", ec="tab:purple", lw=3)
  • reproj vector
  • reproj vector

First option: using the raster as a reference to match, we reproject the vector. We simply have to pass the Raster as an argument to reproject(). See Match-reference functionality for more details.

We can plot the vector in its new projection.

vect_reproj.plot(ax="new", fc="none", ec="tab:purple", lw=3)
reproj vector

Second option: we can pass the georeferencing argument dst_crs to reproject() (an EPSG code can be passed directly as int).

# Reproject in UTM zone 45N.
vect_reproj = vect.reproject(crs=32645)
vect_reproj
Vector(
  ds=             RGIId  ...                                           geometry
       0   RGI60-15.03410  ...  POLYGON ((492957.601 3089314.260, 492981.141 3...
       1   RGI60-15.03411  ...  POLYGON ((494171.838 3088756.723, 494154.107 3...
       2   RGI60-15.03412  ...  POLYGON ((485583.154 3091398.082, 485579.012 3...
       3   RGI60-15.03413  ...  POLYGON ((485023.954 3091501.638, 485023.954 3...
       4   RGI60-15.03414  ...  POLYGON ((484821.599 3093784.463, 484821.599 3...
       ..             ...  ...                                                ...
       81  RGI60-15.10070  ...  POLYGON ((494217.057 3109893.624, 494233.434 3...
       82  RGI60-15.10075  ...  POLYGON ((495932.461 3102375.067, 495918.707 3...
       83  RGI60-15.10077  ...  POLYGON ((496516.356 3105890.705, 496506.863 3...
       84  RGI60-15.10078  ...  POLYGON ((496366.016 3102647.385, 496368.080 3...
       85  RGI60-15.10079  ...  POLYGON ((497125.488 3107721.651, 497091.543 3...

       [86 rows x 23 columns]
  crs=PROJCS["WGS 84 / UTM zone 45N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",87],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32645"]]
  bounds=BoundingBox(left=470889.8334359533, bottom=3080342.6652930346, right=511219.57766997576, top=3113331.061389028))


Total running time of the script: (0 minutes 0.441 seconds)

Gallery generated by Sphinx-Gallery