Skip to content

Instantly share code, notes, and snippets.

@albert-decatur
Last active August 29, 2015 14:19
Show Gist options
  • Save albert-decatur/6b4582c19979bf25fd3f to your computer and use it in GitHub Desktop.
Save albert-decatur/6b4582c19979bf25fd3f to your computer and use it in GitHub Desktop.
see landsat-api scene extents on geojson.io
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
#!/bin/bash
# see Development Seed landsat-api scene extent on geojson.io
# example use: $0 LC81660362014196LGN00
scene=$1
# get scene metadata from devseed landsat-api
curl -s "https://api.developmentseed.org/landsat?search=${scene}" |\
# get corner coords
jq '.results[]|[.lowerLeftCornerLongitude,.lowerLeftCornerLatitude,"",.lowerRightCornerLongitude,.lowerRightCornerLatitude,"",.upperRightCornerLongitude,.upperRightCornerLatitude,"",.upperLeftCornerLongitude,.upperLeftCornerLatitude,"",.lowerLeftCornerLongitude,.lowerLeftCornerLatitude]|@csv' |\
# write as WKT
sed 's:\\::g;s:"::g;s:\(,\{2\}\):|:g;s:,: :g;s:|:,:g;s:^:POLYGON((:g;s:$:)):g' |\
# get geojson from WKT
wellknown |\
# send geojson to geojson.io
geojsonio
@albert-decatur
Copy link
Author

see landsat-api scene extents on geojson.io

sources

  • landsat-api_166036.geojson comes from see-scene.sh, in this gist
  • example landsat scene downloaded with landsat-util

NB

The centroid from landsat-api does not quite agree with the polygon extent generated from API metadata.

Centroid calculated using PostGIS 2.1 w/ GEOS:

db=scratch
for i in *.geojson
do 
    # convert geojson to shp
    shp=$( basename $i .geojson ).shp
    ogr2ogr $shp $i
    # put shp in postgis
    shp2pgsql $shp |\
    psql $db
    # get centroid
    echo "COPY ( SELECT ST_ASTEXT(ST_CENTROID(geom)) AS wkt FROM \"$(basename $shp .shp)\" ) TO STDOUT" |\
    psql $db
done
source wkt
landsat-api POINT(48.5730475299341 34.6062500297351)
landsat-api metadata POINT(48.57291 34.61062)

This table shows distances in meters between these pairs of points.
This was made using a wkt_point_distances.sh, which given a list of WKT points, will output the distances between all non-duplicate pairs.

(Not needed in this case, but handy in general.)

distance_meters first_point second_point meaning
484.939851593 POINT(48.5730475299341 34.6062500297351) POINT(48.57291 34.61062) landsat-api poly centroid to landsat-api metadata centroid

This shows that there is a slight difference between the metadata centroid and the polygon centroid.
However, a quick check in QGIS makes it seem that the see-scene.sh polygon extent and the landsat-util download TIFs almost line up:

full extent
corner zoom - problem!

notes

Example going from 166036 in USGS shp to UTM39N coords:

select st_astext(st_transform(geom,32639)) from ( select st_geomfromtext('POLYGON ((47.336184808754908 34.000773145197954,47.336 34.0008,47.344247813251087 34.031334500589537,47.356167073033426 34.075461184083807,47.723740028525512 35.43626510683103,47.733319293907478 35.471728820126707,47.7439 35.5109,47.744165123892181 35.510861466860142,49.7461 35.2199,49.734986101711954 35.182108728143895,49.315052192874738 33.754181631894454,49.3035 33.7149,47.336184808754908 34.000773145197954))',4326) as geom ) a;

UTM zone can be determined by MTL.
Where can we get these before downloading the scene?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment