Robin's Blog

Converting latitude/longitude co-ordinates to Landsat WRS-2 paths/rows

Updated 6th Jan 2020: This post has been updated to fix the code example and the link to the Landsat shapefile download.

As part of some work I was doing for my PhD, I needed to automatically find what Landsat scene path and row would contain a pixel with a certain latitude/longitude co-ordinate (this was basically so I could automatically download the relevant Landsat image and do some processing on it).

There is an online converter (provided by the USGS) which allows you to input latitude/longitude and get back path/row, but it doesn’t have an API of any kind, and I couldn’t see an easy to way to automatically submit requests to the form, as they’ve (very sensibly really) included methods to stop people doing that.

Summary: So, I wrote my own code to do the conversion. It’s available on Github and depends on GDAL/OGR and Shapely. See the README file on Github for instructions.

Very helpfully, the USGS have made shapefiles of the WRS-2 scene boundaries available to download, so all my code needs to do is load the shapefile, do some ‘Point in Polygon’ operations to find out which polygon(s) the point is in, and extract the ‘path’ and ‘row’ attributes from those polygons. This is fairly easy to do in Python using the wonderful OGR library to access shapefiles, and Shapely to do the ‘Point in Polygon’ operations. When writing this I used the Python Geospatial Development book as a good reference for how to combine GDAL and Shapely in this way.

The code is available at GitHub, with instructions on how to install it (basically, download the shapefiles from the USGS – I’m deliberately not re-distributing them here – and download the Python script), and it is very easy to use:

>>> import get_wrs
>>> conv = ConvertToWRS()
>>> conv.get_wrs(50.14, -1.7)
[{'path': 202, 'row': 25}]
>>> conv.get_wrs(50.14, -1.43)
[{'path': 201, 'row': 25}, {'path': 202, 'row': 25}]

Beware:

  • Instantiating the class may take a while, as it has to load in the shapefile
  • The get_wrs method may return multiple results, as Landsat scenes overlap at all four edges.

Please let me know (via a comment on this post) if you find this code useful.

For reference, the code is shown below, and is well-commented to show how it works:

import ogr
import shapely.geometry
import shapely.wkt

class ConvertToWRS:
    """Class which performs conversion between latitude/longitude co-ordinates
    and Landsat WRS-2 paths and rows.
    Requirements:
    * OGR (in the GDAL suite)
    * Shapely
    * Landsat WRS-2 Path/Row Shapefiles - download from USGS site
     (http://landsat.usgs.gov/tools_wrs-2_shapefile.php), you want wrs2_descending.zip
    Usage:
    1. Create an instance of the class:

        conv = ConvertToWRS()
    (This will take a while to run, as it loads
    the shapefiles in to memory)
    2. Use the get_wrs method to do a conversion:
        print conv.get_wrs(50.14, -1.43)
    For example:
        >>> conv = ConvertToWRS()
        >>> conv.get_wrs(50.14, -1.7)
        [{'path': 202, 'row': 25}]
        >>> conv.get_wrs(50.14, -1.43)
        [{'path': 201, 'row': 25}, {'path': 202, 'row': 25}]
    """
    def __init__(self, shapefile="./wrs2_descending.shp"):
        """Create a new instance of the ConvertToWRS class,
        and load the shapefiles into memory.
        If it can't find the shapefile then specify the path
        using the shapefile keyword - but it should work if the
        shapefile is in the same directory.
        """
        # Open the shapefile
        self.shapefile = ogr.Open(shapefile)
        # Get the only layer within it
        self.layer = self.shapefile.GetLayer(0)

        self.polygons = []

        # For each feature in the layer
        for i in range(self.layer.GetFeatureCount()):
            # Get the feature, and its path and row attributes
            feature = self.layer.GetFeature(i)
            path = feature['PATH']
            row = feature['ROW']

            # Get the geometry into a Shapely-compatible
            # format by converting to Well-known Text (Wkt)
            # and importing that into shapely
            geom = feature.GetGeometryRef()
            shape = shapely.wkt.loads(geom.ExportToWkt())

            # Store the shape and the path/row values
            # in a list so we can search it easily later
            self.polygons.append((shape, path, row))

    def get_wrs(self, lat, lon):
        """Get the Landsat WRS-2 path and row for the given
        latitude and longitude co-ordinates.
        Returns a list of dicts, as some points will be in the
        overlap between two (or more) landsat scene areas:
        [{path: 202, row: 26}, {path:186, row=7}]
        """

        # Create a point with the given latitude
        # and longitude (NB: the arguments are lon, lat
        # not lat, lon)
        pt = shapely.geometry.Point(lon, lat)
        res = []
        # Iterate through every polgon
        for poly in self.polygons:
            # If the point is within the polygon then
            # append the current path/row to the results
            # list
            if pt.within(poly[0]):
                res.append({'path': poly[1], 'row': poly[2]})

        # Return the results list to the user
        return res

Categorised as: Academic, GIS, Programming, Python, Remote Sensing


13 Comments

  1. Beatriz M says:

    Hi,

    Nice post and nice code, thanks! I am also interested on getting Landsat tile bounds from shapefile; the command ‘poly[0].bounds’ would to it in you code.

    Indeed, I tested it. Unfortunately, I noticed that these extracted bounds define a smaller area than the image area. As an example, try to crop a Landsat raster file by passing the extracted WRS coordinates as target boundary; then, display the new file: you will see that you lost some parts of the image. Do you know why is this? Here you have my test:

    $ python
    >>> import get_wrs; conv = get_wrs.ConvertToWRS(); conv.get_wrs(-34.6,146.5)
    [NEW] poly[0].bounds (145.517, -35.510899999999999, 147.928, -33.7149)
    [{‘path’: 92, ‘row’: 84}][{‘path’: 92, ‘row’: 84}]
    >>> exit()
    $ gdalwarp -t_srs EPSG:4326 -s_srs EPSG:32655 -te 145.517 -35.510899999999999 147.928 -33.7149 LE70920842011237ASA00_B6_VCID_1.TIF cropped2WRSbound.TIF
    Creating output file that is 7931P x 5908L.
    Processing input file LE70920842011237ASA00_B6_VCID_1.TIF.
    0…10…20…30…40…50…60…70…80…90…100 – done.

    –> Then display cropped2WRSbound.TIF with qgis/arcgis, you will see what I mean…

    What do you think about that?

  2. Robin Wilson says:

    Glad you found the code useful. I’m not sure why the cropping isn’t working properly. Have you tried doing the cropping manually using something like ArcGIS or QGIS? You should be able to select the right polygon from the shapefile and then use it to crop it – does this give the right result? If it doesn’t, then it suggests there’s something strange with the polygon data themselves – and I’d suggest contacting the USGS.

    If that works ok then it could be that the latitude/longitude you’re running with is one of the points that is within multiple WRS-2 polygons – and in that case you might be using the wrong polygon to get your bounds from. The only other thing that I can suggest is to double and triple check the gdalwarp command: I’ve found in the past that the order of the co-ordinates for cropping isn’t always the logical order, and have run into strange problems. Unfortunately I haven’t really got any other ideas – but if you find out what’s going on then I’d love to hear from you.

  3. Patrick says:

    Hi Robin,
    thanks for the code. But I have a problem and I have no clue whats wrong. I tested the program once and everything went fine. But than something must went wrong because when I do the statement conv = ConvertToWRS() the error comes:

    Traceback (most recent call last):
    File “”, line 1, in
    conv = ConvertToWRS()
    NameError: name ‘ConvertToWRS’ is not defined

    Import get_wrs works fine. I really do not no what happens since I tested bevor and it worked.
    Do you have an idea what I am doing wrong?
    Thanks for your help,
    Patrick

  4. Robin Wilson says:

    Glad you’re finding the code useful Patrick.

    I think the problem might be how you’re importing it: you’ll need to do ‘from get_wrs import ConvertToWRS’ – let me know if that works.

    Robin

  5. Patrick says:

    HI Robin,
    thanks for the fast reply. Sorry, I am a completely newby to python. Now A different message appears. THis is what I have done:

    1) I put the get_wrs.py code into the C:\Python27 folder together with the unzipped wrs2_descending files.
    2) On the Python IDLE shell:
    ‘from get_wrs import ConvertToWRS’
    ‘conv = ConvertToWRS()’

    and this is the message from the Python IDLE:

    Traceback (most recent call last):
    File “”, line 1, in
    conv = ConvertToWRS()
    File “C:\Python27\get_wrs.py”, line 49, in __init__
    self.layer = self.shapefile.GetLayer(0)
    AttributeError: ‘NoneType’ object has no attribute ‘GetLayer’

    Its so funny because when I first tested your code it worked and now it doesnt. The problem remains also when I replace the code with a freshly downloaded version.
    I really wonder whats wrong.
    Cheers,
    Pat

  6. Patrick says:

    Ah, now it works! I specified the absolute path to the shapefile.
    Cheers,
    Pat

  7. Robin Wilson says:

    Glad it’s working now.

    You might find it easier if you don’t put the files in C:\Python27, but run them all from a separate folder. Anyway, glad you’ve got it going.

  8. Karel Jacobs says:

    Hi Robin,

    just needed a simple function to convert from lat/lon to path/row and came up with this:

    https://github.com/22kaj22/nasa_wrs_converter.git

    plan at some point in the future to return multiple results if applicable, but for now it just returns a single tuple.

    Enjoy!

    KJ

  9. Mitch says:

    Hi Robin,

    This is a great tool! In conjunction with the landsat-util python module I’ve been able to automate my landsat search and download from a shape file of points. One problem I encountered is that I need to decide (automatically) which tile to choose when I have overlaps. I encountered a problem where the area I wanted to subset was at the very edge of a scene. I’m kind of new to python so I’m hoping you could help me out. I was thinking about choosing the tile closest to the point I am interested in. Do you know how I can get the centroid location of the tile polygon using your tool? Any help is much appreciated!

    Mitch

  10. Robin Wilson says:

    Hmm, that’s a very good idea Mitch. (One of those things where I wonder why I hadn’t thought of it!). Those centroids should be available relatively easily…I’ll have a play with the code and see how easily I can implement it.

    Cheers,

    Robin

  11. Bang Pham Huu says:

    I’ve added your code to get all the WRS paths which intersect with polygon instead of just point. Thanks for your idea

    https://github.com/bangph/ConvertLatLongToWRS_Landsat_Download/blob/master/README.md

  12. Istivita says:

    Hi, im trying to use your code but it give me error (( ‘NoneType’ object has no attribute ‘GetLayer’ )). I’ve put the absolute path of the shapefile. But actually im not sure about the shapefile, i downloaded it on https://www.usgs.gov/land-resources/nli/landsat/landsat-shapefiles-and-kml-files , and load this >> ogr.Open(‘E:/Spyder/conus_ard_grid.shp’)
    is that correct ?

  13. Robin Wilson says:

    Hi – I’m glad you’ve found my code useful. I hadn’t used it for years, so I just went back and had a look – and I’ve updated the blog post and the Github README to make it actually work now! I think you have the wrong shapefile downloaded – I’ve updated the link to the shapefile in the blog post, and that one should work. Also, you don’t need to open the shapefile in OGR yourself – you can just pass the path to the file to the ConvertToWRS class constructor – like this: ConvertToWRS(shapefile="E:/Spyder/shapefile.shp"). Hope that helps!

Leave a Reply

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