Skip to main content

EarthAI Platform Support Portal

Reprojecting Vector Data

Astraea
  • Updated

What is a CRS or projection?

A coordinate refrence system is a set of mathematical constructs used to translate locations on the three-dimensional surface of the Earth to the two dimensional plane. Knowing what CRS spatial data is expressed in is important because it tells us where on the Earth the data is. Comparisons between spatial data can only be done when both data are in the same CRS.

Sometimes also called a projection, the translation of data from one CRS to another is often called reprojection.

Common Projections

There are several common projections that you will encounter as you work with raster and vector data.

  • WGS84 - Basically the GPS coordinates of the location. If the CRS is not specified, many vector sources can be assumed to be in WGS84. This is also sometimes referred to as EPSG:4326 or CRS84.
  • EPSG:3857 - The most common projection for web mapping services. Also known as web Mercator.
  • Universal transverse mercator projections:
    • Those using the WGS84 model of the Earth:
      • EPSG:326XX for northern hemisphere, for example zone 15N is EPSG:32615
      • EPSG:327XX for southern hemisphere, for example zone 1S is EPSG:32701

Convenience Functions

Let's get started. First, import convenience functions from earthai.geo for reprojecting [Shapely][Shapely] vector geometries in Python.

In [2]:
from earthai.init import *
from earthai.geo import reproject_from_wgs84, reproject_to_wgs84, reproject_on_the_fly, reproject

from shapely.geometry import Point
from shapely.wkt import loads as loads_wkt
import geopandas
import folium
Importing EarthAI libraries.
EarthAI version 1.4.0; RasterFrames version 0.9.0.dev+astraea.1ce1ff3; PySpark version 2.4.4

Creating SparkSession...

  SparkSession is available as `spark`.

Reprojecting a Single Point

Let's define a point in Doha, Qatar. By convention with many of the libraries used in Python and RasterFrames, the longitude is first and latitude is second. Notable exceptions are folium.

In [3]:
doha_wgs84 = Point(51.533333, 25.286667)
crs84 = 'OGC:CRS84'

Reproject to the Qatar national grid. Many countries define one or more CRSs to be used in official maps.

In [4]:
reprojected = reproject_from_wgs84(doha_wgs84, 'epsg:28600')
reprojected.x, reprojected.y
Out[4]:
(231970.459459592, 392639.4879880317)

Reproject it back to longitude and latitude. There is a small numeric error.

In [5]:
reproject(reprojected,
         'epsg:28600',
         crs84
         ).wkt
Out[5]:
'POINT (51.53333299650284 25.28666700325073)'

Reproject a Single Polygon

This is basically equivalent to reprojecting each vertex. Consider a polygon roughly around the border of the US state of Wyoming.

In [6]:
wyoming_wgs84 = loads_wkt('POLYGON((-111.04226573835854 45.03030926505925,-104.05496105085854 45.03030926505925,-104.05496105085854 41.01665169983997,-111.04226573835854 41.01665169983997,-111.04226573835854 45.03030926505925))')

wyoming_wgs84.bounds
Out[6]:
(-111.04226573835854,
 41.01665169983997,
 -104.05496105085854,
 45.03030926505925)

The reproject_on_the_fly function will reproject to the appropriate UTM zone for the geometry, based on its center. It returns both the geometry and the projection that it is in.

In [7]:
g, p = reproject_on_the_fly(wyoming_wgs84)
g.bounds, p
Out[7]:
((-8167.066779450863, 4541035.645358698, 579459.7235281595, 5004106.413363835),
 {'init': 'epsg:32613'})

We reproject back to WGS84.

In [8]:
reproject_to_wgs84(g, p).bounds
Out[8]:
(-111.04226573835854,
 41.016651699839976,
 -104.05496105085854,
 45.03030926505925)

Reprojecting GeoPandas Objects

The official GeoPandas docs covers both setting the CRS for and reprojecting GeoDataFrames and GeoSeries.

Here we'll present a quick example.

In [9]:
gs = geopandas.GeoSeries([Point(-77, 40), Point(-80, 42)])
gs.crs = crs84
gs
Out[9]:
0    POINT (-77.00000 40.00000)
1    POINT (-80.00000 42.00000)
dtype: geometry

Reproject a GeoPandas object using the to_crs method.

In [10]:
gs.to_crs('epsg:32617')
Out[10]:
0    POINT (841487.434 4435426.172)
1    POINT (582818.069 4650259.848)
dtype: geometry

Reprojecting in a Spark DataFrame

In EarthAI Notebook, we have convenient functions for working with columns of vector geometries.

We will start by creating a Spark DataFrame from scratch. The geom column will contain geometry objects.

In [11]:
from pyspark.sql import Row
spark_geo_df = spark.createDataFrame([Row(geom=g) for g in gs])
spark_geo_df
Out[11]:
geom
POINT (-77 40)
POINT (-80 42)

To reproject a column of vector geometries, use the st_reproject function. The arguments are as our reproject function for Shapely geometries: geometry, origin CRS, and destination CRS. However they can each be columns on the Spark DataFrame. For the CRS columns, you can pass a string lit{ open=new}, rf_mk_crs, st_wgs84_crs, or rf_crs.

In [12]:
spark_geo_df.withColumn('geo_3005', st_reproject('geom', 
                                              st_wgs84_crs(),
                                              lit('epsg:3005')))
Out[12]:
geomgeo_3005
POINT (-77 40)POINT (4941964.655833907 874123.66792639...
POINT (-80 42)POINT (4605070.357505446 884976.46880128...

Working with Raster Extents

Most tile columns in RasterFrames have embedded spatial information: the CRS and extent of the tile. The extent can be interepreted as a vector geometry. There are some convenience functions for considering raster extents. The code below shows:

  • querying Earth OnDemand and reading rasters
  • reprojecting the tile extent from the native CRS of the imagery to WGS84
  • converting the WGS84 extent to a GeoPandas GeoDataFrame
  • viewing the result in a web map
In [13]:
raster_df = spark.read.raster(
    earth_ondemand.read_catalog(Point(-112, 33), start_datetime='2019-01-01', end_datetime='2019-05-31',
                                cloud_cover=2, collections='landsat8_l1tp'),
    ['B2']
)
  0%|          | 0/19 [00:00<?, ?it/s]
100%|██████████| 19/19 [00:00<00:00, 91.51it/s]
100%|██████████| 19/19 [00:00<00:00, 54.38it/s]




















rf_geometry lets us extract the native extent as a Polygon object.
rf_crs extracts the CRS of the imagery. It is formatted as a proj string













In [14]:



    


raster_df.select(rf_geometry('B2'), rf_crs('B2'))

    















    
Out[14]:








Showing only top 5 rows


rf_geometry(B2)rf_crs(B2)




POLYGON ((386085 3780135, 386085 3787815...[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs ]


POLYGON ((393765 3780135, 393765 3787815...[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs ]


POLYGON ((401445 3780135, 401445 3787815...[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs ]


POLYGON ((409125 3780135, 409125 3787815...[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs ]


POLYGON ((416805 3780135, 416805 3787815...[+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs ]
























Combine them together to get tile coverage in our desired CRS.













In [15]:



    


extents = raster_df.select(
    st_reproject(rf_geometry('B2'), 
                 rf_crs('B2'), 
                 st_wgs84_crs()
                ).alias('tile_geo_wgs84')
)
extents

    















    
Out[15]:








Showing only top 5 rows


tile_geo_wgs84




POLYGON ((-112.23576968189883 34.1559282...


POLYGON ((-112.15247105836079 34.1567386...


POLYGON ((-112.0691692227752 34.15749259...


POLYGON ((-111.98586440711995 34.1581899...


POLYGON ((-111.90255684342442 34.1588308...
























Finally, we will convert a subset of the DataFrame to GeoPandas and make a map using folium.













In [16]:



    


extents_gdf = geopandas.GeoDataFrame(
    extents.limit(200).toPandas(), 
    geometry='tile_geo_wgs84', 
    crs=crs84
)

folium.Map([33, -112], zoom_start=8).add_child(folium.GeoJson(extents_gdf))

    















    
Out[16]:





















 




          

            
          

        

      


      

      

        
          

        
        
          
  

    
      
Related articles

    
    
  



        
      

      
        

          

            

              

                Comments
              

              
0 comments

              
            


            
    
              
            


            

            

            
Article is closed for comments.