I have to compare geometries anywhere in the world; mainly between 65N and 40S but also in the Arctic latitudes. I need "commercial accuracy" - not military or navigational. For example does a line intersect with a polygon and - if so - so they intersect for 50 km rather than 55 km. I store the lat and lng of  these geometries.   Can anyone please tell me how to select the most suitable cartesian projection (from spatial_ref_sys) to transform the coordinates into before compare any line with any polygon ?  I have tried to use the lng of the centre of the line to select one of the SRIDs between 32601-32660 and 32701-32760. As many of you will know, these are based on the equator every 6 degrees and the 601-660 are for northern hemisphere while the 701-760 are for the southern hemisphere.  To test whether this works, I check whether the ST_Intersects() function detects that a point in Asia at 36.9417, 75.4083 intersects with a polygon of Asia after both are transformed with SRID 32643 whose central_meridian is 75 degrees. Sadly the ST_Intersects() function doesn't detect that the point is right in the middle of Asia, but when I use SRID 32644 whose central_meridian is 81 degrees, ST_Intersects() function detects an intersection between point and polygon.  There are many other examples of my "calculated" SRID being the wrong one. Maybe I'm overcomplicating ? I originally used SRID 4326 anywhere worldwide before testing with ST_Intersects(), and things seem fairly accurate except in high northern latitudes, so I was hoping to minimize inaccuracies by using different cartesian SRIDs at different longitudes.  Perhaps I need to find SRIDs which centre on the latitude of the geometries - not just the longitude ? Maybe I'm overcomplicating ?   Many thanks for your initial welcoming response. I hope that gives good information.
123928
I am trying to customize an Albers and a Hotine Oblique Mercator (HOM) projection to minimuze distortion in the region I am analyzing. The region extends from about 51 to 62 degrees latitude, covering an area about the size of the Ukraine. The region is oriented NW - SE.  I want to make sure I am using correct methods for determining the two projection parameters: lat/long of the projection center, and centerline azimuth.  I am using ArcMap v10.  Here's the procedure I've followed so far:   Created a single polygon that defines the analysis region (by, generally, creating a convex hull around the extent of the watersheds covering the region).  This polygon is the area I am customizing the projection for. Projected the polygon to Geographic/NAD 83. Used Jeff Jenness' Tools for Graphics and Shapes () to determine the polygon's center of mass on the GRS80 spheroid.  The resulting coordinates are what I used for the "projection center" parameter. To determine the centerline azimuth, I first projected the polygon to an azimuthal equidistant projection, specifying the projection center at the coordinates determined in Step 3. Then I drew a polyline (in the azimuthal equidistant projection), snapped to the projection center point, representing the directional trend of the region polygon.  To get the azimuth at the projection's center, I used Jeff Jenness' Tools for Graphics and Shapes to determine the beginning azimuth of the geodesic curve at the central point. For the Albers projection I am using the longitude for the projection center, as determined in step 3. Am also using the awesome spreadsheet created by Bill Huber () to determine where to place the standard parallels to minimize the scale distortion within the polygon region.   If needed to know, I am using ArcMap's version of the HOM that uses a central line defined by one point at the projection's center and its angle of azimuth. ESRI calls this the "Hotine_Oblique_Mercator_Azimuth_Center". In EPSG, I believe this is the Oblique Mercator, Hotine Variant B, EPSG method code 9815.  I am hoping there are some projection experts out there that can tell me if the above procedure, especially Steps 3 and 4, is a correct way to determine the needed projection parameters. Am I on the right track? Is it correct to be determining the center on the spheroid and the angle of the geodesic from the center point (instead of a "2d" geometric center and azimuth)?    I hope the problem description was clear. I'm eagerly looking forward to any answers, tips, discussion, etc.!
111305
I'm fairly new at calculating distances with precision, but I wanted to know what is the best formula for calculating short distances in UTM coordinates.  The distances that I'am calculating are all less than 30 km and all coordinates lie at Central America.   For instance, the pythagorean distance could cause too much bias?   If not, would it still be a respectable choice either way?   If possible, base your answer on some sort manual or academic paper (and make a reference to it).
51535
How would I go about resolving this in C#?  I have a line segment (great circle distance) defined by two Lon-Lat pairs (call them points A and B). A third point somewhere on the Earth sphere.  The line segment closes the shortest distance between A and B. Third Lon-Lat pair is the point C.  How do I calculate the distance between the point C, and the closest point on the AB segment?  EDIT: In fact, I am looking for an algorithm to find the spherical distance between an arbitrary point C and a geodesic segment AB.  What i have for now:  public static float DistanceInKilometres(PointF A, PointF B, PointF C) {     var d13 = Haversine_KM_F(A, C);     var R = 3961.3;     var brng12 = Bearing(A, B);     var brng13 = Bearing(A, C);     var dXt = Math.Asin(Math.Sin(d13 / R) * Math.Sin(brng13 - brng12)) * R;     return (float)dXt; }  public static Double Bearing(PointF coordinate1, PointF coordinate2) {     var latitude1 = coordinate1.Y.ToRadian();     var latitude2 = coordinate2.Y.ToRadian();     var longitudeDifference = (coordinate2.X - coordinate1.X).ToRadian();     var y = Math.Sin(longitudeDifference) * Math.Cos(latitude2);     var x = Math.Cos(latitude1) * Math.Sin(latitude2) -     Math.Sin(latitude1) * Math.Cos(latitude2) * Math.Cos(longitudeDifference);     return (Math.Atan2(y, x).ToDegree() + 360) % 360; }  private static float ToRadian(this float angle) {     return (float)(Math.PI * angle / 180.0); }  private static float ToDegree(this float angle) {   return (float)(Math.PI * angle / 180.0); }
119947
I am unsure of the best approach to this problem. I can think of some ways of doing it which I will list below but I am looking for the best and most efficient way of doing this.  Problem   Given a variable number of points, create a "route" or line between the points. Create a bounding box (of 50km[arbitrary]) around the route. Return all points within this bounding box Points must be returned in order along the route (A->D)   Note: The radius does not have to be circles. A,B,C,D can use rectangular bounding boxes.    Data   Table with points (x,y - technically lat/lon) Input points to create a route Distance from route to include points   We are using PostgreSQL (8.4) and PostGIS (1.3.6) and Python.  This is the solution I came up with. Thoughts? Ideas? Input?  Possible Solution   Create three polygon "tubes" (A->B, B->C, C->D) Filter points to only those in the polygons. Subtract B->C from A->B so the polygons do not overlap (no duplicates) Subtract C->D from B->C so the polygons do not overlap (no duplicates) For each polygon Calculate vector (A->B, B->C, C->D) Find furthest point opposite of vector in the polygon. Calculate distance from furthest point to each site in polygon. Return points in order of distance (smallest -> largest) Join all the points together (A->B points, B->C points, C->D points)   Solution  My approach was naive, after taking some time and looking at the PostGIS methods I came up with this single SQL call. Note its specific to our database but might help someone in the future. Im sure it could be improved as well.      SELECT aerodrome_id     FROM       geo_aerodromecustom_view,        (SELECT ST_Buffer(           ST_Transform(               (SELECT ST_MakeLine(geom)                FROM                    (SELECT geom                    FROM geo_aerodromecustom_view                   WHERE aerodrome_id IN ('CYOW','CYVR')                    ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END                   ) sites               ),               900913           ),           50000) as geo       ) as line     WHERE ST_Intersects(ST_Transform(geom, 900913), line.geo)     ORDER BY ST_Line_Locate_Point(       ST_Transform(           (SELECT ST_MakeLine(geom)             FROM                (SELECT geom                 FROM geo_aerodromecustom_view                WHERE aerodrome_id IN ('CYOW','CYVR')                 ORDER BY CASE aerodrome_id WHEN 'CYOW' THEN 1 WHEN 'CYVR' THEN 2 ELSE 5 END               ) sites           ),           900913       ),       ST_Transform(geom, 900913)     )   I had to project the data because I am using PostGIS 1.3.4 which doesnt support Geography type.  Basically what I am doing is I am  using ST_MakeLine and a query to locate "aerodromes" and return their geometry.  I had to order them (using the CASE directive) so that the line would be connected in the right order.  I then project and buffer this line to create a Polygon that I can then use to see what other aerodromes intersect with the buffered polygon.  Using the unbuffered line (route) and the call ST_Line_Locate_Line I then order the discovered aerodromes as they appear appear along the path.
123781
I've tried this a few different ways and keep getting the same error. I have a shapefile currently in WGS84 UTM Zone 15 which is SRID 32615.  I used the ArcMap multipart to singlepart tool to create a singlepart polygon. I followed the instructions found at  to install postgresql and postgis.  I open pgAdmin III and select the postgis_21_sample table. Plugins > PostGIS Shapefile and DBF Loader 21.  I check connection details and it connects fine.  I add my .shp file. I set Schema to public and Table to postgis_21_sample. I leave geo column alone and set SRID to 3265.  When I click import I get  ============================== Importing with configuration: postgis_21_sample, public, geom, C:\GIS\forest_buffer\postgres_single.shp, mode=c, dump=1, simple=0, geography=0, index=1, shape=1, srid=3265 Shapefile type: Polygon PostGIS type: MULTIPOLYGON[2] Failed SQL begins: "SET CLIENT_ENCODING TO UTF8; SET STANDARD_CONFORMING_STRINGS TO ON; BEGIN; CREATE TABLE "public"."postgis_21_sample" (gid serial, "shape_leng" numeric, "shape_area" numeric, "orig_fid" int4); ALTER TABLE "public"."postgis_21_sample" ADD PRIMARY KEY (gid);" Failed in pgui_exec(): ERROR:  function addgeometrycolumn(unknown, unknown, unknown, unknown, unknown, integer) does not exist LINE 9: SELECT AddGeometryColumn('public','postgis_21_sample','geom'...                ^ HINT:  No function matches the given name and argument types. You might need to add explicit type casts.  Shapefile import failed.   I've also tried projecting this shapefile to WGS84 with SRID 4326 and receive the same error.  There are >100,000 features and I need to buffer, dissolve, negative buffer, and finally buffer back to the original size to fill gaps.  ArcMap and QGIS fail to do it after days of waiting.
91100
I am looking for a projection (or a set of projections) that decently preserve distances so that I can buffer every country in the world and have a similar level of minimal distortion for each.   I then need a projection (or set of projections) to compute the areas of these buffers so that each country in the world's buffer area can be compared with any other country's buffer area.   Eventually I will need to do a bunch of overlays onto these areas to chop up and assign other statistics based on proportion of overlap and such. All of these resulting stats need to be comparable with every other country in the world. At this stage in the game, I am thinking that must be done in an equal-area projection(s) as well.   It is the first time I am operating at a scale smaller than local or regional, and the results will have to be defensible and open to peer-review. I am assuming it is improper to do all of these things in lat/lng, but then when I think about the logic of why, I confuse myself!  I appreciate any info, opinions, reflections on similar experiences, etc anyone has about this!
71807
What is the function for getting SRID from Geometry field.   I have inserted a shapefile into SQL Server table with a Geometry field defined. I know the projection of shapefile is Swreff99_1800 (EPSG:3011).   Now my question is How I get this SRID:3011 by using some OGC method of the geometry field?
44771
I have a table of Roads in PostGIS which geometry is "geometry(MultiLineString,4326)".  Its attribute are      gid   osm_id  name    ref type    oneway  bridge  maxspeed    geom   Now i want to know the closest Road from latitude and longitude.  Also I want to set some distance from that Line/Road like 20m both side and if the point is within that distance then it should give the name of that line else not.
4203
While importing my shapefile data to PostGIS, I did not select the proper Projection.  How do I now change the SRID of the data, without transforming the Coordinates?
113120
I have a road network and want to test connectivity between the features.  I am having trouble with two specific features, the rest of the dataset is functioning just fine.  When I get those 2 links' start &amp; end locations using ST_StartPoint and ST_EndPoint, the result is that they are slightly different.  However, when I convert the result to text by wrapping it like so:  ST_AsText(ST_StartPoint(geom)) and ST_AsText(ST_EndPoint(geom))  then the results are exactly the same.  Is it possible they are actually off by a tiny amount like .000000001 and the ST_AsText loses a decimal point?
44745
I have been attempting to use the Spatial Analyst Curvature tool with mixed results.  I originally tried to use it with a raster that was projected using a 1983 Teale Albers projection and the results were black, indicating the lowest values.  I then transformed the raster into UTM projection and ended up with grid-like lines on the raster which have values that are obviously not accurate.  Does anyone know why these lines have appeared or what I'm doing wrong.  Thanks.
80614
I'm trying to find the points within a certain distance, in meters, from another point.  I'm using the function ST_Buffer, but I don't get how to express the distance in meters.  I'm trying with something like this:  select 1 as cartodb_id, ST_Buffer(   ST_Transform(     ST_GeomFromText(       POINT(-58.38145 -34.60368)'       , 4326     )     , 3857   )   , 500 ) as the_geom_webmercator   I have the coordinates from google maps (in wgs84, 4326) then I transform them to webmercator (3857) and then I pass 500 as parameter, and I get something that looks like the desired area, but after comparing it with what google maps says it falls short for a non neglectable distance.  Is this the right way to achieve it?   note: you can try the sql query in cartodb  --  From this question:  I came with the following code (I use SRID 31997 to calculate the buffer, and then go back to webmercator)  select 1 as cartodb_id, ST_Transform(   ST_Buffer(     ST_Transform(       ST_GeomFromText(         'POINT(-58.427185 -34.624592)'         , 4326       )       , 31997     )     , 2000   ), 3857 ) as the_geom_webmercator   Now if falls short for no more than 20 meters. I thought mercator (EPSG 3857) was in metric units.
52632
There are several algorithms to calculate distance between two latitude/longitude points.     Hubeny     But there are not so many algorithms (actually I can't finally find it) which can calculate latitude/longitude of point B from another latitude/longitude point A with bearing and distance.      I'm glad to know if there are any algorithms which is lesser accurate, but faster than Vincenty's formulae.
33403
I want to do something relatively simply but it seems to be very difficult to do, though I am very new to using R for spatial data so apologies if I've missed something obvious.  I have a list containing 161 shapefiles of bird distributions in the world, they are each in the format of a SpatialPolygonsDataFrame. I also have world climate data in a raster stack format for 8 variables. Both sets of files are in the same projection (lat and long WGS84). I want to create new columns for the climate data with 1's where the species polygon overlaps with the climate data and 0 where it doesn't. I tried various functions to do this, quite a few initial options didn't work, the best option appeared to be:  pnt.in.poly(XY coordinates of climate data points, XY coordinates of species polygon)   I changed my raster of climate data into points to get the XY coordinates. I unified the polygons for each species shapefile and then used:  species@polygons[[1]]@Polygons[[1]]@cords   to get the coordinates for my species polygons as a matrix. However, I discovered that:  species@polygons[[1]]@Polygons[[1]]@coords    wasn't actually getting all the coordinates of the polygons only some of them.  I tried using fortify from the ggplot package, and it gives me a dataframe with coordinates but some of them are for holes in the polygons. I ran that with the pnt.in.poly function excluding the holes but it took an hour to run and when I looked at the output maps after it matched for some of the range but not all of it- there were holes in the distribution where there shouldn't be holes.   I also tried this method that I found somewhere else:  lonlat &lt;- lapply(slot(species, "polygons"), function(x) lapply(slot(x,"Polygons"), function(y) slot(y, "coords")))    which gives me a list with 436 items- but I'm still not sure how to get the coordinates into a matrix form that will work in pnt.in.poly.  I'm sure this should be easier to do! Does anyone know a better way?
107072
I'm having difficulty wrapping my mind around meridian convergence, the horizontal angular distance between true north and grid north.   Say I want to determine if a feature lies to the east or west of the central meridian (CM) by using true north and grid north designations on a topo map.   I realize grid north is parallel to the CM, while true north points to the poles, and magnetic north is the direction a compass points and constantly changing.   Looking at the figure below, I see grid north points to the east of true north. So does this fact alone mean everything on the grid is to the east of the true north, and therefore the CM? Does this simple fact mean that features in this grid lie to the east of the CM, or am I missing something.   Also, I realize I can just look at the tick marks and and see that this this grid in particular covers -80°22’ to -80°15, and the CM for zone 17R lies at 81°. But had I not known that and wanted to use the north arrows. Please explain.
105586
I've been looking over the PostGIS "API" of functions, and I notice that most of them take two elements to compare. For example, the ST_Distance function takes two geometry/geography elements to find the distance.   There is no function to do something like: "Given a geometry G, give me the closest geometry GClosest in Table T where G.id &lt;> GClosest.id"  I realize I could write a PL/PgSQL function to iterate over the table and call ST_Distance on each element, but I'm hoping there's a better, more efficient, solution.
54545
I'm using PostGIS 1.5.2, with geometries in SRID:900913. I need create a circles using a list of points as centre, with a radius of 600 kilometers. I'm using this query:  INSERT INTO circles (geom) (    SELECT ST_Buffer(point, 600000, 'quad_segs=8')     FROM points );   But circles created have not 600 kilometers of radius (radius is near this length, but not exactly).  Are there other methods to create circles in PostGIS?  NOTE: Information represented is from Spain. Correct projection is 4326, but client use Google rasters, so I am storing data in 900913 to avoid reprojections and increase performance.