10

I have just started working in PostGIS. I have a table in PostgreSQL that stores Latitude,Longitude. I want to create a rectangle of 500*500 foot around every point. Therefore, used the below given script.

INSERT INTO GridData(latitude,longitude,geom)

SELECT 
  latitude,longitude,
  ST_Buffer(ST_SetSRID(ST_MakePoint(longitude,latitude),4326),152.40)
from public.PointData

I am using Geometry data type in GridData and PointData table.The above script has been written to create rectangular shape in GridData table. What is the correct way to do this?

M.Sharma
  • 313
  • 2
  • 11
  • Welcome to GIS SE. Thank you for taking the Tour. A rectangle in geographic coordinates isn't going to look rectangular. You're using ST_Buffer, which doesn't generate rectangles. Your ST_MakePoint syntax appears to have the lon/lat flipped (unless you're mapping in Antarctica within 30nm of Amundsen–Scott Station). Please [Edit] the question to specify whether you're using geometry or geography datatype and where your data exists, and choose one question to ask. – Vince Jul 12 '17 at 11:59
  • Hi Vince, I have updated my question. Please help. – M.Sharma Jul 12 '17 at 12:12
  • You did not specify where your data is. – Vince Jul 12 '17 at 12:19
  • @Vince : My data is in Postgresql PointData table that will be populated with Latitude-Longitude from some location of USA. I want that a rectangle box should be created around every point and then copy those boxes in GridData table. – M.Sharma Jul 12 '17 at 12:27
  • You are curating a question. Do not place thanks in it. – Vince Jul 12 '17 at 12:35
  • Last edit: What does the current command produce? To me it looks like a constant circle for all rows. Hint: Use longitude and latitude in your expression, then decide if you want to use ST_Extent or ST_Envelope as per this question: https://gis.stackexchange.com/questions/14948/bounding-box-for-postgis-table – Vince Jul 12 '17 at 12:43
  • Current command produces rectangle.But not sure about there size that I required. However, first I will go through the link you gave me and check the result. I do not know how to attach screenshot , otherwise I show you result. – M.Sharma Jul 12 '17 at 12:58
  • @Vince: you are correct my current command does not product rectangle. I was wrong they are just point. I go through the link https://gis.stackexchange.com/questions/14948/bounding-box-for-postgis-table but this did not solve my problem. – M.Sharma Jul 13 '17 at 00:04
  • I was able to use ST_Expand for a similar problem. – fqhv Nov 20 '19 at 01:58

3 Answers3

7

The ST_Buffer function will buffer the specified distance in the units of the geometry's spatial reference system. In your case, you are effectively asking to buffer a distance of 152.40 degrees, which is not really what you want to do.

Also, buffering around a point will create a circle and not a rectangle. To get your rectangle you can instead use the ST_MakeEnvelope function which requires you to specify your minimum and maximum x and y extents which you can determine from your point coordinates.

If you want to work in units of feet or meters you will need to project your points to a projected coordinate system that uses feet or meters. This coordinate system will depend on the extent of your data; if it covers the contiguous USA you could use the US National Atlas EPSG:2163, or if it is more regional you should look into the US State Plane projections and UTM Zones. Of course you also have to consider the accuracy you are hoping to achieve.

The following SQL example constructs a point geometry in WGS84 before transforming to the US National Atlas projection and determining the minimum and maximum x and y extents. Since the point is going to be the centre of the rectangle these extents are calculated using half of the intended length/width. Additionally, since the US National Atlas projection uses meters I've converted the 250 feet to the corresponding distance in meters.

SELECT 
ST_MakeEnvelope(
    (ST_X(ST_Transform(ST_SetSRID(ST_MakePoint(longitude,latitude),4326),2163))-(250/3.28)), --Min X
    (ST_Y(ST_Transform(ST_SetSRID(ST_MakePoint(longitude,latitude),4326),2163))-(250/3.28)), -- Min Y
    (ST_X(ST_Transform(ST_SetSRID(ST_MakePoint(longitude,latitude),4326),2163))+(250/3.28)), -- Max X
    (ST_Y(ST_Transform(ST_SetSRID(ST_MakePoint(longitude,latitude),4326),2163))+(250/3.28)),  --Max Y
    2163
    )

There is certainly a much more efficient way of constructing the SQL but hopefully this illustrates the points above.

Ali
  • 4,055
  • 1
  • 17
  • 29
  • I use this for my testing and will tell you the result soon. I have a question that you use EPSG:2163 when use ST_Envelope, but what will happen if coordinates are also from outside USA (For example canada ,UK ,Germany or hong kong). Then do I need to change this. if yes, then doesn't there any common projection? – M.Sharma Jul 13 '17 at 13:48
  • 1
    @M.Sharma If your coordinates are outside of the projection's intended area (USA in this case) they will be subject to greater distortion. You certainly don't want to use 2163 for the UK. The topic of selecting a coordinate system is a complex one with many considerations - the following page is worth a read: https://gis.stackexchange.com/questions/2769/what-strategies-criteria-or-rules-to-use-for-selecting-coordinate-systems. For greater accuracy you can project each point to its corresponding UTM zone before making the envelope. – Ali Jul 13 '17 at 14:08
5

For generating a square for each points we can use the combo ST_ENVELOPE(ST_BUFFER()):

For example if we want a square of 500 units (so we need a 250 units buffer):

SELECT ST_Envelope(ST_Buffer(MyGeom,250))

Which is not really efficient since we need to compute some buffers, but it is short and readable.

obchardon
  • 1,724
  • 2
  • 13
  • 21
1

For my needs (not an exact 500 foot square, 10% is close enough, so adjust the division to suit you), this is another way to create the square envelope around a LINESTRING that runs from bottom-left to top-right.

Note that at this latitude, there are approximately 100000 metres per degree (111.3 km at the equator).

WITH myconstants (Longitude, Latitude) as (
   values (28.0, -26.0)
)
select ST_Envelope(ST_GeomFromText(CONCAT('LINESTRING(', 
                                          Longitude - 250.0/3.28/100000.0, ' ', 
                                          Latitude - 250.0/3.28/100000.0, ',', 
                                          Longitude + 250.0/3.28/100000.0, ' ', 
                                          Latitude + 250.0/3.28/100000.0, ')'), 4326)  ) FROM myconstants
PeterS
  • 668
  • 5
  • 16