2

Given a point, latitude and longitude, how would one find the coordinates of a rectangle centered on that point with a known width and height in meters? I tried using the Geodesic stuff from geographiclib in python to compute the rectangle coordinates, but the resulting rectangle has the wrong aspect ratio. The general idea here is to compute a rectangle around a given point and then use that rectangle to clip LANDSAT 8 imagery to a final image with a known aspect ratio.

Sorry if this is a dumb question, I've just started doing work with GIS recently.

[EDIT]

I tried modifying my code to do something like what the first answer below does, but this doesn't quite work for me, I run into the same problem I was having before, when I clip the landsat tiff I get an image with the wrong aspect ratio. The command I'm using to clip the geotiff looks something like this gdalwarp -of gtiff -t_srs EPSG:4326 -te <bounds from the polygon that I calculated> input.tif output.tif. The width and height I used should give me an image with an aspect ratio of 1.6 but instead I get an image with an aspect ratio of about 1.9 or so.

  • 2
    What about to convert lat, lon to plan coordinates (e.g. using web mercator projection or to the projection that will be used to display images) and define there the rectangle. – Zoltan Feb 27 '16 at 6:01
3

Here's some example code using pyproj. Given a point in lat lon, it calculates new lat lon points given a distance in meters and an azimuth. The azimuth comes from the aspect ratio of the rectangle.

from math import sqrt,atan,pi
import pyproj
geod = pyproj.Geod(ellps='WGS84')

width = 10000. # m
height = 20000. # m
rect_diag = sqrt( width**2 + height**2 )

center_lon = -78.6389
center_lat = 35.7806

azimuth1 = atan(width/height)
azimuth2 = atan(-width/height)
azimuth3 = atan(width/height)+pi # first point + 180 degrees
azimuth4 = atan(-width/height)+pi # second point + 180 degrees

pt1_lon, pt1_lat, _ = geod.fwd(center_lon, center_lat, azimuth1*180/pi, rect_diag)
pt2_lon, pt2_lat, _ = geod.fwd(center_lon, center_lat, azimuth2*180/pi, rect_diag)
pt3_lon, pt3_lat, _ = geod.fwd(center_lon, center_lat, azimuth3*180/pi, rect_diag)
pt4_lon, pt4_lat, _ = geod.fwd(center_lon, center_lat, azimuth4*180/pi, rect_diag)

wkt_point = 'POINT (%.6f %.6f)' % (center_lon, center_lat)
wkt_poly = 'POLYGON (( %.6f %.6f, %.6f %.6f, %.6f %.6f, %.6f %.6f, %.6f %.6f ))' % (pt1_lon, pt1_lat, pt2_lon, pt2_lat, pt3_lon, pt3_lat, pt4_lon, pt4_lat, pt1_lon, pt1_lat)

The documentation for pyproj.Geod can be found here.

Below is a screenshot of the point (yellow) and rectangle (green) in QGIS: enter image description here

  • This doesn't quite work for me, I run into the same problem I was having before, when clip the landsat tiff I get an image with the wrong aspect ratio. The command I'm using to clip the geotiff looks something like this: gdalwarp -of gtiff -t_srs EPSG:4326 -te <bounds from the polygon i obtained> input.tif output.tif. The width and height I used should give me an image with an aspect ratio of 1.6 but instead i get an image with an aspect ratio of about 1.9. – Paul Wicks Feb 27 '16 at 19:05
  • Ok, hopefully we can chase this down. So the code I provided gives the same undesired result? Does the clipping polygon itself have the correct aspect ratio? Does the landsat tile completely contain the polygon? I can see how the result could have the wrong aspect ratio if they do not completely overlap. You could drop the tiff into QGIS and use the QuickWKT plugin to make the polygon. Maybe visualizing them at the same time might give insight into the problem. – user55937 Feb 27 '16 at 23:33
  • Ooops, messed up on that last comment, here's my actual reply: >> So the code I provided gives the same undesired result? Yes >> Does the clipping polygon itself have the correct aspect ratio? No, it has the same incorrect aspect ratio I think. >> Does the landsat tile completely contain the polygon? Yes, I checked that already. – Paul Wicks Feb 28 '16 at 0:18
  • Glad we ruled out the polygon generation step. How are you measuring aspect ratio? I wonder if the aspect ratio of the result is measured in pixels, i.e., width/height of the pixels in the clipped image. And if the pixel resolution is tied to degrees, then it won't seem to have the same aspect ratio as it would in linear units. I.e., squares aren't squares when draped over spheres kind of thing. Can you try to convert it to UTM and see what the aspect ratio is with a linear unit, or post the image so I can? – user55937 Feb 28 '16 at 0:57
  • Maybe plug your clipping coordinates into here: rcn.montana.edu/resources/converter.aspx to get the utm zone. Convert your landsat to UTM using gdal_translate, and try to clip the image in a linear CRS to see if the image actually has a different aspect ratio than you expect. – user55937 Feb 28 '16 at 1:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.