texts: What are the set of projections that can be accurately georeferenced using polynomial 2?

texts:
Background I'm starting with a map from the internet (Encyclopaedia Britannica) and georeferencing it using QGIS 2.0.1, I get really accurate georeference solution using polynomial 2 transformation: very small residuals, most less than 5 pixels. (This is surprising because the previous time I tried this, I got pretty poor results: Improving Georeferencing results? though now I am using more points (31) and was quite careful with them.)

The georeference window:

(note the very slight residuals, the biggest of which is the small line emanating from Beijing). This produces the following overlay, which you can tell is quite good.

Goal I'd like to identify the projection and its parameters used by the cartographer.

Questions Can the fact that a polynomial 2 fits really well be used to narrow down the list of projections that the map uses? Or do many/most projections georeference well with polynomial 2 when applied over such a large geographic extent and with so many control points?

(I can disable most of the control points and as long as I include at least one in the northern regions, I can get the same very good georeference. And of course polynomial 3 works very well also.)

(I've also tried to fit most of the pseudocylindrical projections mentioned in The Proj Manual (1996, revised 2003) to see if one can be made to match the control points as well as the polynomial 2 georeference fit, but the latter is much better than any of those fitted projections.)

The GCPs I have put the GCPs at http://pastebin.com/pWti0W8k. The map x/y are in EPSG:4326 (I used the Google Physical OpenLayers map as the reference but then reprojected those EPSG:3857 to EPSG:4326 so I could get lon/lat).


texts: How to change legend orientation in QGIS?

texts:
I created a legend in the print composer and the default layer orientation is

How can i change the legend elements orientation like in ArcGIS?

I seeked in the Item properties but didn't found any hint.


texts: Percent of data pixels in a raster clipped to a buffer in PostGIS?

texts:
I have a raster created by unioning a set of rasters that are clipped to a buffer. I'd like to determine whether this "unioned raster" has a sufficient number of data (vs. no-data) pixels. I'm wanting to ensure my understanding of the unioned raster and how I'd compute this is correct.

The way I am doing it now is to get the width (st_width) and height (st_height) of the unioned raster, and compute the area of an ellipse of that (since the raster is not square in dimension). This in theory gives me the total number of possible pixels in the buffer that this raster represents. Then, I get the actual total count of data-pixels by doing a sum on the st_valuecount from that unioned raster.

Some code snippets per the above. Note, that I have many tables I'm doing this with, so while I could wrap all this up in a function, I can't since I can't dynamically substitute table names in (to my knowledge), but I imagine there's a more efficient way to get that unioned raster and use it for all these.

Get the width (and similar for height):

WITH myconstants (buffer) as (
    VALUES (st_transform(st_buffer(latlng_point_with_srid(latitude, longitude, upgis_utmzone_wgs84_from_latlng(latitude, longitude)), buffer_radius), 4326))
)
select st_width(st_union(st_clip(rast, buffer)))
    FROM TABLE-NAME-GOES-HERE, myconstants where st_intersects(rast, buffer);

Then do π x (0.5 x width) x(0.5 x height) for the area/total possible pixels.

Get the count of pixels with data:

WITH myconstants (buffer) as (
    VALUES (st_transform(st_buffer(latlng_point_with_srid(latitude, longitude, upgis_utmzone_wgs84_from_latlng(latitude, longitude)), buffer_radius), 4326))
)
select sum((pvc).count) from (select st_valuecount(st_union(st_clip(rast, buffer))) as pvc
    FROM TABLE-NAME-GOES-HERE, myconstants where st_intersects(rast, buffer)) as x;

And then I could do something like count/area as a measure of what percentage of pixels have data.

As an example with concrete numbers, I have a raster that reports a width, height of 311, 210. It also reports a count of data pixels of 50534. So:

raster total area: 311 x 210 = 65310
approx possible total pixels: π(156 x 105) = 51459
percentage of data pixels: 50534 / 51449 = 98.2%

Is this correct, or what may I not be understanding or doing incorrectly, etc.?


texts: Connect to an ESRI Personal Geodatabase (*.mdb)

texts:
My environment:

Windows 7 Professional:
Service Pack 1
64-Bit

QGIS:
QGIS-Version: 2.8.2-Wien
32-Bit

Microsoft Office 2010 Professional:
32-bit

I want to connect an ESRI Personal Geodatabase (*.mdb) with the "Add Vectorlayer" feature in QGIS

The problem is with the Access Jet driver.

My Questions:

Has anybody successfully connected to an ESRI Personal Geodatabase in QGIS 2.8.2?


texts: How to work with Windows SMB connections in QGIS on Linux/Debian?

texts:
I'm having trouble accessing our Windows shared drives in QGIS while working on Linux Mint.

For some reason I don't remember this being a problem prior to version 3.4, but here is the issue:

Our shared team drive sits on a Windows server, which we access via a mapped drive in Windows. This works well in QGIS on Windows, as the standard 'browse' dialogue boxes show mapped drives, etc.

However, the browse dialogues in QGIS on Linux - now using Mint v19 - don't show any detail:

Can anyone think of a way to access drives that have been added via SMB and available at the desktop level, but not accessible to QGIS?


texts: Exporting Georeferenced Models from Meshlab?

texts:
I am struggling to understand Meshlab's and 3D mesh file's relationships to georeferencing. I am coming from a GIS background, so I believe my assumptions on how this should work are preventing me from succeeding with 3D software.

I've found plenty of good tutorials for georeferencing in Meshlab online, but I haven't found one that discusses exporting the mesh with embedded coordinates to move to another software package with the benefit of the georeferencing.

I have an OBJ file of an archaeological site that I generated using photogrammetry software that is not georeferenced. I am using a text file of coordinates and the georeferencing dialog in Meshlab to georeference the mesh using the "pick a point on moving" method. This generally works great. And when I save the project and reopen it the mesh displays as georeferenced.

However, my expectation based on the way GIS works and based on how I figure scaled models should be exchangeable between software suites is that there is some way to write the georeferencing information to the mesh itself (rather than saved in the Meshlab project file).

I've exported the georeferenced mesh as OBJ and PLY but when importing into a new Meshlab project the model is no longer georeferenced. I've tried simplifying the coordinates (removing the leading digits on the UTM coordinates to emulate a smaller scale coordinate system) that I read can help if the coordinates are larger than the software is anticipating. But this hasn't made a difference.

There are numerous reasons I would benefit from the mesh itself being georeferenced. For example loading multiple georeferenced models into Meshlab at a later date to view their relationships, upload a georeferenced model to SketchFab so that is is properly scaled, and to create a stand alone scaled deliverable for a client or colleague that may use different 3D software.

Is there a separate command from "Export" to write this information to the mesh?


texts: Error 000666, does not have OID with Locate Features Along Route in ModelBuilder

texts:
I'm trying to take a bunch of point data and determine what points of surveyed lines are closest to each point.  To this end I have put together a model that loops through the line features and in each loop it selects the survey points that intersect the line, creates Thiessen polygons from those points, then uses Locate Features Along Routes to generate a table showing what is covered where.  The problem is that the Locate Features along Route tool consistently (but not universally, for some reason) spits out an error (#000666) saying that the input features do not have OIDs.

I've tried the Make Query Table method, but that doesn't work either, since the Thiessen polygons are wiped each time through the loop the Query Table loses the definitions of what fields to take, resulting in a different error, since the Make Query Table method will spit out a table, rather than a shapefile, if you don't check the option to include the shape in the input.

Does anyone know of another workaround to give those polygons OIDs?  Or of another way to get a similar result without using the clearly buggy Locate Features Along Routes?


texts: CartoDB or Leaflet.js get coordinates of point clicked?

texts: Using CartoDB or Leaflet.js is it possible to retrieve the coordinates of the point on the map clicked when opening an infowindow?


texts: Errors installing rgdal on LINUX system?

texts:
I am trying to install rgdal on a Linux cluster, but am running into some issues. I have successfully installed gdal and proj4, but when I try to build rgdal, I receive the message that "proj_api.h not found in standard or given locations". proj_api.h is located at /bin/proj4/include/proj_api.h. I have added the path to $PATH, however I still can't compile it.

How can I specify the path so that I can install rgdal?


texts: Finding closest outside point to point inside polygon in PostGIS?

texts:
There is the famous ST_ClosestPoint() return the closest point inside polygon to some outside point.

I'm looking the inverse problem: I have a location inside polygon, and I want to get the closest location outside the polygon.

I started thinking that first step should be to find the closest point on the boundary of the polygon. Although I don't know how to achieve this, I'm not sure how to continue even from this stage (considering I have the closest boundary point to my inside point).

EDIT:
What about comparing results from ST_Line_Locate_Point() running on all lines of the polygon and choosing the minimum: This gives me the closest point on the border line

WITH border_closest_points AS (
        SELECT ST_Line_Locate_Point(border_lines.geom, ST_GeomFromText('POINT(lon lat)', 4326)) AS fraction,
        ST_Line_Interpolate_Point(border_lines.geom, ST_Line_Locate_Point(border_lines.geom, ST_GeomFromText('POINT(lon lat)', 4326))) AS closest_point_on_line
        FROM border_lines
)
SELECT border_closest_points.closest_point_on_line AS geom, ST_Distance(ST_GeomFromText('POINT(lon lat)', 4326), border_closest_points.closest_point_on_line) AS distance
FROM border_closest_points
ORDER BY distance ASC LIMIT 1;


texts: How to select a given feature in a GeoJSON file using Python?

texts:
Scenario:

I have a large (120+ MB) .geojson file that represents a bunch of LineStrings features.

Here is an example feature selected using json.dumps(geoFile['features'][0]):

{
    "type": "Feature",
    "id": 0,
    "properties": {
        "FID": 0,
        "prop1": 1,
        "prop2": "thing2",
        "prop3": "thing3"
    },
    "geometry": {
        "type": "LineString",
        "coordinates": [
            [
                -99.491525,
                35.306851,
                0
            ],
            [
                -99.491485,
                35.306705,
                0
            ],
            [
                -99.491451,
                35.306581,
                0
            ],
            [
                -99.491417,
                35.306456,
                0
            ],
            [
                -99.491375,
                35.306321,
                0
            ]
        ]
    }
}

Goal:

Select a feature based on the property value in prop1.
Render the selected feature on a geopandas map (or other python map) along with some georeferenced points from dataframes being used elsewhere in the app.

Questions:

What is the Python syntax for selecting the feature from the GeoJSON file based on properties/prop?

Obviously I know how to select based on feature index ([0] above), but I can't find examples of selecting by prop1.
Tried: json.dumps(geoFile['features'][0]['properties']['prop1'] == 1) but this just returns true


texts: Directional pattern in normalized LiDAR point cloud

texts:
I have used the lidR package (https://github.com/Jean-Romain/lidR) to produce a normalized point cloud. I then filter out the points tagged with the overlap bit using:

opt_filter(ctg) <- "-drop_overlap"
catalog_retile(ctg)

The resulting point cloud appears to have a directional pattern, where points are present in somewhat uniform strips. It seems the unfiltered point cloud has a similar pattern but much it is much less pronounced. Is this a normal feature of lidar data?

Here is a 3d representation of the normalized and filtered point cloud with the pattern I have described:


texts: QGIS 3.4 server on Windows 10: enable logging without fcgid?

texts:
I'm running QGIS 3.4 server on Windows 10.  The docs here say that you need to add these lines to httpd.conf to enable QGIS server logging:

  FcgidInitialEnv QGIS_SERVER_LOG_FILE /var/log/qgis/qgisserver.log
  FcgidInitialEnv QGIS_SERVER_LOG_LEVEL 0

(Those lines are also referenced in several other: here and here for example)

But those only work if mod_fcgid is installed.  Apparently you have to jump through some moving-target hoops to install mod_fcgid on Windows.

Is there a way to enable QGIS server logging without mod_fcgid?


texts: Inaccurate Transformations - OSGB 1936 / British National Grid to WGS84 with Proj.Net

texts:
I am new to all this GIS work, so if I have done something silly or just plain wrong please just let me know.

I m trying to covert from coordinates from OSGB 1936 / British National Grid to WGS84 but I am seeing different results to other tools.

My example coordinates X,Y is; 362895, 155602

My conversions appear to be nearly correct however they are always a little out. Here is a link showing the conversion on epsg.io. If you switch the result to decimal you will see the expected results is;

Long=-2.5335813°, Lat = 51.2983258°

https://epsg.io/transform#s_srs=27700&t_srs=4326&x=362895.0000000&y=155602.0000000

I have also done the same conversion on https://mygeodata.cloud/cs2cs/ with the same results as epsg.io.

However in Proj.Net my result is;

-2.53225739401083, 51.2985133921956

Here is a code example to show my issues;

class Program
    {

        public const string wktOsgb1936 = "PROJCS[\"OSGB 1936 / British National Grid\"," +
        "GEOGCS[\"OSGB 1936\"," +
        "DATUM[\"OSGB_1936\"," +
        "SPHEROID[\"Airy 1830\",6377563.396,299.3249646,AUTHORITY[\"EPSG\",\"7001\"]]," +
        "AUTHORITY[\"EPSG\",\"6277\"]]," +
        "PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]]," +
        "UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]]," +
        "AUTHORITY[\"EPSG\",\"4277\"]]," +
        "UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]," +
        "PROJECTION[\"Transverse_Mercator\"]," +
        "PARAMETER[\"latitude_of_origin\",49]," +
        "PARAMETER[\"central_meridian\",-2]," +
        "PARAMETER[\"scale_factor\",0.9996012717]," +
        "PARAMETER[\"false_easting\",400000]," +
        "PARAMETER[\"false_northing\",-100000]," +
        "AUTHORITY[\"EPSG\",\"27700\"]," +
        "AXIS[\"Easting\",EAST]," +
        "AXIS[\"Northing\",NORTH]]";

        static void Main(string[] args)
        {
            // The orginal 27700 OSG co-ords
            var x = 362895;
            var y = 155602;

            var res = TransformSpatial(wktOsgb1936, x, y);

            // the long/lat which I am expecting as confirmed on other systems
            var shouldx = -2.5335813;
            var shouldy = 51.2983258;

            // A small calc to see the difference
            var difx = res[0] - shouldx;
            var dify = res[1] - shouldy;

            Console.WriteLine("From    : {0}, {1}", x, y);
            Console.WriteLine("conv    : {0}, {1}", res[0], res[1]);
            Console.WriteLine("Should  : {0}, {1}", shouldx, shouldy);
            Console.WriteLine("Diff    : {0}, {1}", difx, dify);

            Console.ReadLine();
        }

        public static double[] TransformSpatial(string wktFrom, double x, double y)
        {
            var cf = new CoordinateSystemFactory();
            var f = new CoordinateTransformationFactory();
            var sysFrom = cf.CreateFromWkt(wktFrom);
            var transformTo = f.CreateFromCoordinateSystems(sysFrom, ProjNet.CoordinateSystems.GeographicCoordinateSystem.WGS84);
            var ret = transformTo.MathTransform.Transform( x, y );
            return new double[] {
                ret.x, ret.y
            };
        }

    }

Please help and let me know what I am doing wrong, this is doing my head in!!

Is there a problem in my "From" mappings, my "To" mapping, is there another step I am missing or is there a bug?

Thanks in advance


texts: Automaticallly select and replace all occurences of a font face

texts:
I have around 50 mxd file that is done some times earlier. The text font used in labeling was not licensed. Now I can Change the text in layout element. But it has also some text in map view too. But arcpy.mapping.ListLayoutElements can not find the text document.

Is there any solution in arcpy?


texts: Populating attribute table with attributes from table using PyQGIS

texts:
I am trying to populate an attribute table of a feature with attributes from another feature. I have been able to segregate the input data into fields named "From", "Eastings" and "Northings". I find this is not such a difficult task to carry out as this code helped. The problem now is the input data has only four(4) rows but when I run my script, it multiplies the number of rows and I end up with twenty rows. I am quite sure I need to change something in the code, i am not sure where as I am relatively new to python.
here is my code:
pr.addAttributes([QgsField('FROM',  QVariant.String)])
vpr.addAttributes([QgsField('Eastings',  QVariant.Double)])
vpr.addAttributes([QgsField('Northings',QVariant.Double)])
vlayer.updateFields()

vlayer.startEditing()
for things in vlayer.getFeatures():
    for i in new_attribs:
        things['FROM'] = i[0]
        things['Eastings'] = i[1]
        things['Northings']=i[2]
        things.setAttributes(i[0])
        vpr.addFeatures([things])

        vlayer.updateFeature(things)

    vlayer.commitChanges()

This is the result I get:

How do I get only the four(4) rows I need?


texts: Where to find ArcGIS (Desktop / Server) application icons list (not map symbols)

texts: I tried this on Google but couldn't find anything! Where can I find list of ESRI Applications Icons (not Map Symbols) representing the ESRI concepts like (Layer. Shapefile, Geodatabase, ....)?


texts: Automatically select small polygons contained within currently selected large polygons in QGIS?

texts:
I am using QGIS (2.4.0) and I am trying to identify regions in Australia that are close to working mines. For this purpose I have 2 layers:

Layer 1 - 535 polygons defining the various regions across Australia
Layer 2 -  Points indicating the location of operating mines

I have used "select by location" to identify the polygons in layer 1 that contain a mine.

However my problem is that there are some small polygons in layer 1 (based around towns) that don't actually contain a mine but are within a large polygon that does. I would like to find a way to select these polygons automatically.

Ultimately I am looking for a function that does something like: "Select polygons in Layer 1 entirely contained within currently selected polygons (in layer 1)".

The image below shows one example of what I want to do. The large polygon is already selected using select by location because it contains a mine. I then want to be able to automatically select the smaller polygon contained within it (and do this across the whole layer).


texts:
I have set up a basic QGIS project combining various WMS services and want to publish it to my company's intranet for everyone to use but I am a relative beginner at GIS and all the guides I have found are far too jargony or for previous versions so don't seem to apply.

I have tried the QGIS2Web plugin which was easy to use but didn't give me the popups, and the Lizmap plugin which I thought would be the answer but the user guide skims over uploading so having set it all up I can't work out how to get it from a QGIS file to something that is viewable on the web (presumably html?). By doing some Google searches I found a guide that suggested I needed Lizmap Web Client so I downloaded that but that was talking about putting it in Apache2 (whatever that is) rather than it just being something you install which confused me.

Can anyone direct me to a more beginners guide for creating a webmap on either of these platform? Specifically, is it possible to get QGIS2Web to show popups or is that just not going to happen? Or how do I get my Lizmap QGIS files on to my FTP in a format that is visible online?

Alternatively does anyone have any suggestions for a different approach to take?

Apologies if this seems a bit generic but I am a complete beginner so can't really be more specific than this!

I am using QGIS 3.8.2-Zanzibar, QGIS2Web 3.7.0, and Lizmap 3.1.3.


texts: I want to display some markers on Openlayers map mobile version. I know that standard way of adding markers on Openlayers does not work with Openlayers mobile library. The example of Openlayers site didn't help me add a marker dynamically. Can somebody show me the correct way of implementation?


texts: Trying to write a script that re-projects vector datasets in a folder

texts:
I am new to Python and am struggling with my code specifically with my loops.The project is supposed to abide by these rules:

Must re-project shapefile vector datasets in the folder to match the target dataset's projection.
Must append "_projected" to the end of each projected dataset name. For example: CityBoundaries_projected.shp.
Must skip projecting any datasets that are already in the target projection.
Must report a geoprocessing message of which datasets were projected. In this message, the dataset names can be separated by spaces. In the message, do not include datasets that were skipped because they were already in the target projection.
Must not contain any hard-coded values such as dataset names, path names, or projection names.
Must be made available as a script tool that can be easily run from ArcToolbox by someone with no knowledge of scripting.

I keep getting a this error:

OSError: "C/Users/jessi/Desktop/CPSUGeog485/Lesson2" does not exist

This is my code so far:
import arcpy
import os

#establish workspace
folderWorkspace = "C/Users/jessi/Desktop/CPSUGeog485/Lesson2"
targetFc = "C:/Users/jessi/Desktop/CPSUGeog485/Lesson2/Countylines.shp"
featureClass = ["CountyLines", "Ferries", "PopulatedPlaces", "CityBoundaries", "StateRoutes"]

#establish spatial reference of selected feature class
targetDesc= arcpy.Describe(targetFc)
targetSr = targetDesc.SpatialReference
targetSrName = targetSr.Name

arcpy.env.workspace = folderWorkspace
fcList =  arcpy.ListFeatureClasses()

# List current feature classes in folder
for fcCurrent in featureClass:
    fcCurrentDescribe = arcpy.Describe(folderWorkspace)
    fcCurrentSr = fcCurrentDescribe.SpatialReference
    fcCurrentSrName = fcCurrentSr.Name

for fcCurrentSrName in fcCurrent:
    if fcCurrentSrName != targetSr:
        print ("Error Matching Spacial reference")
else:
    print ("Spatial reference matching succesful")

for fcCurrentSrName in fcCurrent:
    if fcCurrentSrName == targetSrName:
        continue

else:
    # Output feature class and projected name
    fcOut = os.path.splitext(fcCurrent)[0]
    arcpy.Project_management(fcCurrent, fcOut, targetSr)
    print (fcOut +"_projected.shp")```


texts: Attribute displays shape that is not supposed to be there

texts: I have an issue, that a few shapes (in this case Russia and Australia) has a line that spreads accross the globe. I cannot select the line individually and when i select Russia, the line by Russia is automatically also selected. In the attribute table, it just shows up as one shape. How do I solve this?


texts: TMS tileOrigin on Google Maps ImageMapType


texts: What (where) are elements that comprise WGS60 and FSHR1960?

texts:
Considering two 1960 geodetic ellipsoids (from this list), I've found the following parameters (from here):

WGS60:   World Geodetic System of 1960: 6378165.0, 298.3
FSHR1960: Fischer – 1960 (Mercury): 6378166.0, 298.3

This is a new subject for me. With my current understanding, I believe that along with the two ellipsoid parameters, each of these should also include a list of locations or landmarks, and the coordinates of those locations in ellipsoidal coordinates. For each, the collection of the ellipsoid parameters together with the list of locations and their coordinates is called a datum.

If I am correct so far, are the lists available somewhere?

So far I have found many references to these, but not an actual complete listing of either datum.


texts:
I'm having trouble understanding the concept of a Vector tile. My understanding of Vector is point, line, polygon. So buildings, roads, areas etc.

My understanding of Raster is a grid of pixels each with a value to describe it. So the actual imagery would be raster.

I don't understand how a client can request a vector tile and the "imagery" (which I'm used to being raster - .tiff,.png,.jpeg etc) be presented as a vector.

How are these raster images being "vectorized"?

How can points, lines and polygons create something as complex as a raster image?


texts: Linking numerical data to county geometric data using CartoDB?

texts:
I have a CSV spreadsheet that contains some county names in Michigan and numerical data. It looks like this:

Allegan, 5, 4.51
Alpena, 1, 3.44
Antrim 2, 8.64
Barry 1, 1.69

Etc.
I am trying to make it into a choropleth. I went through the tutorials on CartoDB which allowed me to make a choropleth out of data they gave me that already had geometric data for all US counties (including Michigan counties!). My CSV does not have the geometric data. I'm very new to this. I'm not sure how to proceed. I'm assuming either I

extract the Michigan county data from the US county data and link it to my dataset (no idea how to do this)
Hand draw all 83 Michigan country maps myself (seems redundant but willing to do it if this is the only way)
Something else?


texts: gdal_translate of Tiff to gpkg with only one band

texts:
I try to use gdal_translate in my windows command prompt to convert a single band .tif raster to a single band .gpkg.
gdal_translate -of GPKG example.tif example.gpkg

I always get a .gpkg with four bands.
I tried all possible combinations of the following commands with no succes:
-band
-mask
-expand

gdalinfo of example.tif
Driver: GTiff/GeoTIFF
Files: example.tif
       example.tif.aux.xml
Size is 5334, 5280
"
"
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
"
"
Band 1 Block=5334x1 Type=Byte, ColorInterp=Gray
  Min=10.000 Max=113.000
  Minimum=10.000, Maximum=113.000, Mean=50.602, StdDev=14.032
  NoData Value=255
  Metadata:
    STATISTICS_MAXIMUM=113
    STATISTICS_MEAN=50.60229149027
    STATISTICS_MINIMUM=10
    STATISTICS_STDDEV=14.032043952498

gdalinfo of example.gpkg:
Driver: GPKG/GeoPackage
Files: example.gpkg
Size is 5334, 5280
"
"
Metadata:
  AREA_OR_POINT=Area
  IDENTIFIER=example
  ZOOM_LEVEL=5
Image Structure Metadata:
  INTERLEAVE=PIXEL
"
"
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha


texts: arcpy: get list of adjacent lines per line

texts:
I am working with polylines with hundreds of thousands of features, and would like to write a script that returns the OIDs of all lines adjacent (intersecting with) each line. I have been experimenting with the Geometry object and the shp.touches() method, but it seems as though a cursor with a nested cursor on the same feature is the only way to use it, and this takes much too long.
The main goal is to perform an "Eliminate" task on a polyline, where short lines are merged to adjacent longer lines that share some unique ID. Collecting information in a dictionary has always been my technique of choice because of its speed, but dealing with geometries complicates this significantly. Unfortunately, Eliminate only works on polygons, and is quite slow itself.
I have already tried to record shape geometries in a dictionary (dict[row.OBJECTID] = [row.Shape_Length, row.Shape, row.uniqueID]) and attempted a nested iteration, but the geometry objects don't seem to be interacting in the same way they do within a cursor iteration (shp.touches() doesn't work). Has anyone come across any techniques I could use for this problem?


texts:
Is it possible, that my created WPS (in Geoserver) returns an image or a pdf (general: returns a datatype, that isn't a supported format by Geoserver).

Is it sufficient to expand the ProcessParameterIO class?

Can anybody give me a general overview, how to do that?


texts: Creating map book via ArcGIS data driven pages with a dynamic table (graphic table element) falling outside of layout extent

texts: I am using ArcGIS to create a map book that is data driven and has a dynamic table (graphic table element) associated to each page. The data driven page is indexed to follow a line for its extent. Each extent has a different number of points that are represented through the dynamic table. The question I may ask is, is there a way I can export the maps to pdf to include parts of the table that might fall outside of the layout extent (rollover) on a separate page?


texts: FME Desktop - Combining two weighted fields into one

texts:
I have two feature classes that I have spatially joined using the spatial relator. There are two city fields in this new feature. For the sake of the example, we can refer to them as Cities A and Cities B. Some values are null in both fields, however the Cities A is more accurate.
How do I combine Cities A and Cities B into a new field with priority given to Cities A, assuming it is not null. i.e. if Cities A is null and Cities B has a value, we use Cities B. If Cities A and B both have values for a feature, then only Cities A would be used.


texts:
last summer I made the relief of a cave, and at this time I returned with Therion software. This survey is in 3d and I can save it in the formats .Lox, .3d, vtk.
I would like to import it in grass or qgis, software that I know better than others. I can not find the command that allows me to do so. Can anyone suggest me something?


texts: Exporting OpenStreetMap building blocks boundaries to polygon shapefile

texts:
I am using QGIS software to draw building blocks boundaries that appear on OpenStreetMap as shown in the following image

but as you can see it can not be accurate.

Is there a way to export the building blocks boundaries from OpenStreetMap to a polygon shapefile??

I already checked the website Map data and also tried OSM downloader plugin in QGIS but it does not have such required data.


texts: Is it possible to reference data using mapped network drives instead of UNC path in arcgis server services?

texts:
Has anyone managed to create a service using datasets referenced via a mapped network drive letter on windows?
ESRI support claims it is not "supported", but I want to make sure it is actually not possible or just "not supported" as in "not recommended".

the permissions are properly set up for my SOC account, since the services created using UNC path to reference datasets are working fine. I guess it all comes down to making the ArcGISSOC account understand those letters.  I think I tried logging on as the ArcGISSOC and mapping the drives manually, but it doesn't seem to work.

Any ideas?

Thanks


texts: Is there a tool that is similar to the ArcGIS repair geometry tool in QGIS?

texts: I need a tool in QGIS that has the same functionality as the repair geometry tool (native of ArcGIS).


texts: Issue with using arcpy.GetParameter and arcpy.GetParameterAsText to write layer source to text file

texts:
I am trying to get a Map Document and (a new) text file as user defined parameters, in order to write to the text file a list of layer sources. The text file is then suppose to open once complete.

When I hard code in a test mxd and txt inputs the code works, however when I use arcpy.GetParameter/arcpy.GetParameterAsText and set up the code within a Toolbox script, the code doesn't work. It will not create the text file, or even populate a test text file (if hard coded).

I have also tried: outputTxtFile = arcpy.GetParameter(1) with no success.

I am wondering if I am using arcpy.GetParameter/arcpy.GetParameterAsText correctly?

Note: in the Toolbox script properties, I have the first parameter set to 'ArcGIS Map Document' and the second set to 'Text file'.

# Import system modules
import arcpy, os

try:
    #Set input parameters
    mxd = arcpy.GetParameter(0)
    outputTxtFile = arcpy.GetParameterAsText(1) #r"D:\10_GIS\Projects\Maps\Test2.txt"

    #Open txt file
    outputTxtFileWrite = open(outputTxtFile, "w")

    lyrList = []

    for lyr in arcpy.mapping.ListLayers(mxd):
        if lyr.supports("DATASOURCE"):
            lyrSource = lyr.dataSource
            lyrList.append(lyrSource)

    lyrList.sort()
    for i in lyrList:
        outputTxtFileWrite.write(i + "\n")
        arcpy.AddMessage("Writing to file")

    outputTxtFileWrite.close()
    os.startfile(outputTxtFile)

except Exception as e:
    print arcpy.GetMessages()
    arcpy.AddMessage(e)


texts: Raster import from GRASS to R - handling nan/NA

texts:
I am generating a regression model using several environmental covariates. Before importing the GIS layers to GRASS I had them as ESRI raster format. I exported them to ASCII file, and from ASCII file I imported them to GRASS. I did this a couple of momths ago, I believe, because I was having trouble finding a different way of importing the layers directly from ESRI raster.

In GRASS I did several operations and corrections of my raster files. After, I used the package spgrass6 to import the raster files to R, directly. I created a SpatialGridDataFraME

In my preliminary data exploration I see a lot of nan, and I wonder whether I did any mistake exporting the files to ASCII files, or whether is these nan's are normal, since I have many empty cells in my grid (cells outside my study region, but that are within the box, or the grid defined by the region in GRASS).

I want to get rid of NA values, and so I tried:

region.pt <- as.(region, 'SpatialPointsDataFrame')

But it gave me an error:

  cannot allocate vector of size 2.3 Gb

My SpatialGridDataFrame is 1.9 GB!!!

Should I clean all these nan values BEFORE importing them into GRASS, or before importing them into R?

Or is it OK to oblige the spatialGridDataFrame to become SpatialPointsDataFrame?

What will happen with nan's?

How can I get only the pixels within my study area?


texts: Converting csv to kml/geojson using ogr2ogr or equivalent command line tools?

texts:
I've just unsuccessfully tried to use ogr2ogr to convert CSV into KML like so:

ogr2ogr -f "KML" output.kml input.csv

The csv has "latitude" and "longitude" columns.

The attributes are all loaded but if you inspect the KML's text, the Coordinate elements are not being added so its not working.

I feel like this will work if I can inform gdal to interpret certain columns ("latitude", "longitude") as coordinates.

How might I do this?


texts: Trigger selection event for Arcgis Operations Dashboard with js

texts:
I have several dashboards which i want to display within the same webpage.
each dashboard is loaded within an iframe and switching between them is done with basic js and css.
they are all manipulated with the same dropdown category selector using the same layer to help users select the data they want to see.

what i want to do is carry the selection over between all the dashboards, meaning that when the user want to see the indicators for city A he could select it in dashboard 1, and dashboards 2,3,4 and 5 will show the data for city A too.

I tried triggering click and focus events for the  elements, as well as removign and adding the is-active and highlight classes for them like suggested here.

I also tried setting the values in the input and the value-container divs.

I realise that  shouldnt be activated by click events so I tried using the advice from here.

It seems like i'm missing an event that could be triggered, but i can't find any more advice about how make this selection.

Clarification

I know how to access the elements within the <ifram>, i'm looking for a way to create a slection.
here's a codepen with an example dashboard within an iframe and the way to access the elements within the dashboard document.
https://codepen.io/bogind/pen/mYOpKB


texts: Write output of QGIS Batch Processing to temporary layers

texts: Is it possible to write the outputs of QGIS Batch Processing to temporary/memory layers?


texts: Raster to GeoTIFF file conversion problem

texts:
An extension to this question, so when I export the raster as a .tif file, the resultant GeoTIFF file is not similar to the plot in R below. Am I missing a step in the following code? See the images below for comparison. The respective files can be downloaded from here.
d <- read.csv("sample_data.csv")
coordinates(d) <- ~X+Y
data<-data.frame(d$X,d$Y,d$VIC_SC_04)
spdf <- SpatialPixelsDataFrame(d, tolerance = 0.016129, d@data)
# Create the raster
 r <- raster(spdf, layer=3)
# Set the projection to WGS-84
crs(r) = "+proj=longlat +datum=WGS84"

# Export the raster
rstr_WGS84 <- writeRaster(rstr,filename = "~/Sample_raster.tif", overwrite=TRUE)


texts: How to efficiently extract specific metadata for many shapefiles?

texts:
I have a directory with 50 subdirectories containing a total of 300 shapefiles which I would like to query for their feature count and output to a text file.

Thus far I am able to get the feature count for all shapefiles in 1 subdirectory using ogrinfo:

for %f in (*.shp) do ogrinfo -al -so %f

This outputs a lot of additional metadata and tags which need to be stripped out manually.

Introducing the following to this script results in the feature counts all equaling 1:

-dialect sqlite -sql "SELECT COUNT(DISTINCT('ID'))"

I'd love to learn a command line solution to efficiently parse the ogrinfo output. I recognize that this could require some awk/grep type solutions which I can learn to do, but I am working with a standard Windows 10 install of OSGEO /GDAL.


texts: Invoke Google Earth address search from a hyperlinked command in a table cell

texts:
I would like to be able to launch a Google Earth search in Windows by clicking a table cell that contains a street address (plus town and US state abbreviation). For example, I have a table in ArcMap that has a field populated with Google Earth web search urls like this:

https://earth.google.com/web/search/123 Main St, Anytown, NY

The corresponding feature layer can be set so that when you click a feature with the info tool, you can then click the search attribute and launch a Google Earth web search (assuming Chrome is set as the default browser).

But what I would like to find is a way to populate that field instead with a command such that when I click the cell in the table itself, it invokes a Google Earth desktop search on the given address. This because I want to see the search location in the context of kml data I have loaded in Earth on my computer.
I get the impression from this Google Earth Engine site that this might be possible, but there I'm straying outside my competencies.


texts: Can I extract value from a polygon that contains countries to Grid polygons that I have created?


texts: Creating polygons layer from point layer with conveks buffer

texts:
I have a point layer of >750000 points representing dwellings with district heating.

Now, I want to create polygons for every n of points with a distance of n meters to the next point.

I have a polygon layer from a municipality, but the polygons are wrong in some places.

The points outside the dotted-polygon, should have a polygon of themselves. Note this is only a very small fraction of the data.

Also there are polygons where there are no points intersecting.

So, my question is:
Is there somehow a function or something else, which can create the polygons for every eg. >=5 points with a distance og eg. 10 meters to the next point?

I have tried working it out with the convex hull function, but it doesn't do the trick, because I have to select the points by hand, it that process will take ages for the whole dataset.

I'm using QGIS 2.6


texts: QGIS Georeferencer - Image file imports too small

texts:
I am trying to georeference a largish aerial image (50MB) to some GIS data I have and import it as a background.  However the image imports extremely small - so small I can't see it!

I am using QGIS and the Georeferencer plugin.  I select about 7 control points randomly across the image and I am using the thin plate spline with nearest neighbour resampling.  The error is 0.00768.  I have tried several transformation settings and set the target CRS to the project CRS, but get similar results. The control points appear in QGIS in the right positions.

I would appreciate any assistance.

ADDED - A screen shot of the settings used.


texts: trying to update empty string field with a string

texts:
I need to update a column of values such that empty fields are changed to the string "OFF_ROAD."  I've looked up a plethora of examples, and this is what I'm trying to get to work, but when I run it, nothing happens to the input file-

import arcpy

from arcpy import env

env.overwriteOutput = True

arcpy.env.workspace = "F:/EMILY/TEST/GET_TYPE_TEST/"

infc = "SpJoinscript.shp"

cur = arcpy.UpdateCursor(infc)

for row in cur:

if row.TYPE_1 == "":
    row.TYPE_1 = "OFF_ROAD"
cur.updateRow(row)
del row

del cur

print "done mothaf"

what am I overlooking here???!!!


texts: Where to download qgis 3.0 for Mac to avoid the 3.2 bug?

texts:
So the editing function is great and all but I had to get a new laptop and installed this 3.2 which had the nasty bug of only using half the screen

To circumvent this I'd need to download 3.0 and it's not available on their web page. Is there like an alternate repository for Mac? On kingchaos page they only have up until 2.18


texts: Displacing colliding point symbols in QGIS

texts:
I have a layer that includes points for electric poles (Circles) and transformers (triangles). Transformers can be on poles or on their own. I would like to displace the transformers only if they are on top of a pole, keeping their exact location when they are on their own so they follow the electric wires. I have tried a conditional in the displacement expression box:
if( intersects(  $geometry ,  get_feature(  @layer_name , "Type" , 'Pole' )),1,0)

but it is not working.

What I don't want

What I want ONLY when there is a Pole (circle) present.


texts:
I am looking for help to determine a coordinate system used in 1981 in south east BC, Canada. The coordinate system is retrieve from some old logs, and I need to plot and find  the point locations. An example of the data is:

23, 707.3 N
14, 491.6 W

I have no idea what this coordinate system is. I need to figure out how to plot this data into UTM 11.


texts: Export shapefile to CSV files by a grid I drew in a image

texts: The shapefile is in orange and I need to export them into (1.csv and 2.csv etc...) corresponding to each grid I drew in the map. Is there any way this can be done in a single process?


texts: How to convert an IGeometryCollection into a Polyline

texts:
I am using the CostPathAsPolyline method which returns the Polyline as an IGeometryCollection.  I now need to get this converted back into a polyline.  But it seems to be more complex than just a cast.

Can someone please help me?  I cannot find anything online.


texts: export an aerial image to .SID file and a .Tiff file

texts:
I am using arcmap 10.4.1 ( also have PRO 2.0 is this is an valid option)

I have an aerial image that was broken down to small grids and need to export each grid to a .sid file and a .tiff file. I cannot find a tool that supports both of these file types.

Can this be done in arcmap or do I need another software?


texts: Overlay layers on different displays from the same coordinate system

texts: I have vector layers from digitizing in Google Earth and vector layers from digitizing of a pdf document. All the layers were projected in the same UTM zone 30N system (EPSG code 32630). However, layers obtained from Google Earth do not overlap with layers obtained from document. What can I do to fix the problem with these layers overlapping. I work with QGIS 3.10.


texts: How to query all branches of a specific bank in a spcific city

texts:
I would like to know how to query for an example all the amenities that are named Deutsche bank in a city named Berlin in Germany.
Below is my attempt.
Furthermore, please let me know the meaning of the following sign
<  >

Code:
{{geocodeArea:Germany}}->.searchArea;
(
  nwr(area.searchArea)["amenity"="bank"];
  >;
  nwr(area.searchArea)["brand"="Deutsche Bank"];
  >;
  nwr(area.searchArea)["atm"="yes"];
  >;
  nwr(area.searchArea)["addr:city"="Potsdam"];

);
out;


texts: Print esri Map doesn't show feature layer

texts:
I Have a map with some feature layer. I make a print function to print map to pdf and it run correctly.
But when i set filter to feature layer, the layer doesn't show on the pdf print. Please help for my problem.

This the sample code :

 app.flroad.setDefinitionExpression("Year=" + "'" + right(periode[0], 4) + "'" + " and month=" + "'" + left(periode[0], 2) + "'");
    app.flbound.setDefinitionExpression("Year=" + "'" + right(periode[1], 4) + "'" + " and month=" + "'" + left(periode[1], 2) + "'");


texts:
I want to copy and paste (or export) data from a query (or attribute table) into excel.
I've tried copying and pasting direct from a query results (in DB Manager) or loaded the query result as a layer, and done the same from the attribute table.
But in excel, if there is a linefeed or carriage return character (or both) in any cell of the query results, this gets split across multiple rows in excel.
Is there a way to control or avoid this behaviour?


texts: autocad-map of real estates with own coordinate system and information to be converted to use it in QGIS

texts:
I have a question regarding CAD and QGIS, but unfortunately I don't have a clue of CAD.

Here it goes: Real estates data has been mapped. There are several reference points surrounding the area. So there is an own coordinate system being used. The coordinates of the corners of each polygon are known. Each polygon holds different information regarding size, owner etc.

The task is to have this finished mapping converted to be used in QGIS. The person wants me to have it correctly transferred into UTM-PCS. The outcome is an editable map of the real estates and each polygon is given it's equivalent attribut table.

Any hints???


texts: Simplify geometry of multiple features without losing shape

texts: I am simplifying administrative boundaries using a high simplify tolerance (which is my intention so as to get a minimalist appearance).  However, after simplification, the borders of the regions no longer fit together.  What would be the most effective way to solve this problem?


texts:
I got two sets of points.

The first set has a unique ID
The second set references to this ID (one-to-many).

So I merged the two sets. Now every point in the joined table has the X and Y coordinates from the related first set.

Name;ID;X;Y;X_2;Y_2;Distance

Every X and Y is unique but X_2 and Y_2 repeats from a given set.

Now I would like to calculate the distance between (X,Y) and (X_2,Y_2) in the new column for each row.

I tried distance() but I'm way too new in Python to get it working.

I also tried some tools like Distance Matrix, but I got results for all distances (or to the number defined) but not only to the distances in one row. Or between Points sharing the same ID.


texts: I have a table with county polygon boundaries (named "counties"). I need to build a query that identifies a particular county as well as any adjoining counties that share the same border. Have tried breaking this down but I'm missing something. Any ideas?


texts:
Given a degraded image of a Google maps tile, for example:

Is there a service or technique I can use to find the map coordinates of this tile?


texts: Why is transparancy of GRASS rasters missing when imported into QGIS?

texts:
I have a raster map in QGIS imported from a GRASS mapset with the QGIS GRASS plugin. The problem is that I have null() cells in the GRASS map that appear "colored" in QGIS. How can I set the transparency?

I have tried in 'Proprieties' but it seems that it doesn't work. I would like to know a simple way to solve my problem.


texts:
I am trying to georeference rasters from Vietnam in QGIS which have printed coordinates in VN2000.

I am using the following custom CRS settings for the province.

  +proj=tmerc +lat_0=0 +lon_0=108.25 +k=0.9999 +x_0=0 +y_0=500000 +ellps=WGS84 +towgs84=-191.90441429,-39.30318279,-111.45032835,0.00928836,-0.01975479,0.00427372,0.252906278 +units=m +no_defs

The rasters appear in the incorrect location when georeferenced.

This has been solved. The definition of the false easting above is incorrect. The definition should read +x_0=500000 +y_0=0


texts: Use image to apply unique filterDate range for each pixel in an imagecollection

texts:
I have 2 images of timestamps (in system:time_start type format, ms since unix epoch) corresponding to a time series image collection; each pixel in the two images represents a date range, first image is the start date, second is the stop date.

Is there a way to reduce the image collection where each pixel in the image collection is filtered to its own unique time window?


texts: Osm2pgsql failed due to error when I try to import data into PostgreSQL

texts:
I am using the following link to install openstreetmap server on my ubuntu 14.04 server. I have followed their tutorial line by line and installed most of the dependencies and everything. When I try to import the North America's map data into my PostgreSQL database it is throwing an error. I am using this command to load all the data

osm2pgsql --slim  -d gis -C 16000 --number-processes 1 /usr/local/share/maps/planet/north-america-latest.osm.pbf

but it is throwing an error as you see on image.

Osm2pgsql failed due to ERROR: Error: Connection to database failed: FATAL: role 'azureuser' does not exist.

azureuser is the username of my terminal and when configuring my postgresSQL database I created a user 'openstreet' and database 'gis'. If you open the link you can also see the command to import data under the section Importing data into the database but that is not working for me. Could anyone help me with this error? I am not able to move forward from here


texts: Directly generate XYZ tiles from GeoJSON?

texts:
I would like to directly generate XYZ (=raster) tiles from a ~20 GB large GeoJSON file.

This is way too much for QGIS Tiler.
As another trick I've generated vector MBTiles from the GeoJSON, but sadly mbutil only produces garbage tile PNGs.
GeoServer randomly omits objects at higher zoomlevels, but I need to see ALL building at every zoomlevel.

Do you know any other solution for directly generating XYZ tiles?


texts: Why is select features button not working in ArcGIS 10.3?

texts:
Every time i opened a mxd file (in ArcGis 10.3) the default of most of the layersis in "Not Selectable". I search a way to change the default in the table of content ,so all the layers will be always in "selectable" default like the layer "gvul":

I try to change it in selection options and make sure that interactive selection method is "create a new selection" but i didn't find a solution.


texts:
I have tile server running based on mapnik, mod_tile, and renderd. It was succesfully generate and serving tile using mapnik xml configuration which defined on renderd.conf. Something like this :

[default]
URI=/school_tiles/
TILEDIR=/var/lib/mod_tile
XML=/home/brain/school/map.xml
HOST=localhost
TILESIZE=256

Now the scenario level up. I must renew the map based on data which queried by user. Example : if i hit http://xyz.com/?q=school&date=2014-01-01 then mapnik should regenerate map based on that query string and server response back with new tiles.

Is there any possibilities to change mapnik xml in renderd.conf on the fly ?, or what the best technique to solve this scenario?

Thanks in advance.


texts: Upload specific or all tables from Geopackage to PostGIS database

texts:
I would like to upload a specific table (i.e. U2018_CLC2012_V2020_20u1 (Multi Polygon)) from a Geopackage file to a PostGIS database from the command line. I tried:
ogr2ogr -f "PostgreSQL" "PG:database" \
    -t_srs EPSG:3035 -overwrite -lco GEOMETRY_NAME=geom -nln "corine.corine12_vec" \
    U2018_CLC2012_V2020_20u1.gpkg

However, I just realized, the Geopackage contains multiple tables and only (the wrong) one (i.e. U2018_CLC2012_V2020_20u1_FR_MYT (Multi Polygon)) is uploaded.
Q: How can I upload a specific table or all tables?
Output of ogrinfo U2018_CLC2012_V2020_20u1.gpkg:
INFO: Open of `U2018_CLC2012_V2020_20u1.gpkg'
      using driver `GPKG' successful.
1: U2018_CLC2012_V2020_20u1 (Multi Polygon)
2: U2018_CLC2012_V2020_20u1_FR_REU (Multi Polygon)
3: U2018_CLC2012_V2020_20u1_FR_GLP (Multi Polygon)
4: U2018_CLC2012_V2020_20u1_FR_GUF (Multi Polygon)
5: U2018_CLC2012_V2020_20u1_FR_MTQ (Multi Polygon)
6: U2018_CLC2012_V2020_20u1_FR_MYT (Multi Polygon)


texts: How to make Habitat Suitability Map from rasters QGIS

texts:
We're trying to make a HSI (Habitat Suitability Index) model map for a species. The thing is we have a bunch of vector layers corresponding to parts of an entire area and we have to assign a score to each one of them. So in each layer, I created a new attribute and assigned a score in each layer between 0 and 3.

Then I rasterized all those layers, but how can I join all rasters to make a map with all scores by colors?

I can only do it in vector layers with style-categorize- color by attribute


texts: Show WMS information in popup on Leaflet map?

texts:
I need to view WMS information through a popup in Leaflet map. I tried with the code below, but the popup doesn't come out with information, but automatically downloads a file.
map = L.map(document.getElementById('mapDIV'), {
center: [-19.000, 35.09],
zoom: 5.6,
});

var ProvinciasWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Provincias_Moz',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
});

var DistritosWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Distritos',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
});

var PostosAdministraticosWMS = L.tileLayer.wms('http://localhost:8080/geoserver/Mocambique/wms?',{
    layers:'Mocambique:Moçambique_Distritos_e_Postos',
  version: '1.1.0',
   transparent:'true',
   format: 'image/png',
   attribution: "myattribuition"
});

map.addEventListener('click', onMapClick);
    popup = new L.Popup({maxWidth: 1000});
    function onMapClick(e) {
    var latlngStr = '(' + e.latlng.lat.toFixed(3) + ', ' +         e.latlng.lng.toFixed(3) + ')';
    var BBOX =         map.getBounds()._southWest.lng+","+map.getBounds()._southWest.lat+","+map.getBounds()._northEast.lng+","+map.getBounds()._northEast.lat;
    var WIDTH= map.getSize().x;
    var HEIGHT = map.getSize().y;
    var X = map.layerPointToContainerPoint(e.layerPoint).x;
    var Y = map.layerPointToContainerPoint(e.layerPoint).y;
    var URL = 'http://localhost:8080/geoserver/Mocambique/wms?SERVICE=WMS&VERSION=1.1.0&REQUEST=GetFeatureInfo&LAYERS='+'Mocambique:Provincias_Moz,Mocambique:Distritos,Mocambique:Moçambique_Distritos_e_Postos'+'&QUERY_LAYERS='+'Mocambique:Provincias_Moz,Mocambique:Distritos,Mocambique:Moçambique_Distritos_e_Postos'+'&STYLES=&BBOX='+BBOX+'&FEATURE_COUNT=5&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&X='+X+'&Y='+Y;
    popup.setLatLng(e.latlng);
    popup.setContent("<iframe src='"+URL+"' width='400' height='200'         frameborder='0'></iframe>");
    map.openPopup(popup);
    }


texts:
I need to setup 2 SQL triggers in CartoDB that will update my inland and coastal tables each time my obs table is updated.  I have never setup a SQL trigger but I have looked at the documentation and it appears as though I need to have a function stored somewhere that the trigger activates.

So far my first trigger looks like this:

CREATE TRIGGER update_inland
AFTER INSERT ON obs
FOR EACH ROW
EXECUTE PROCEDURE update_inland()

But I get an error message that "function update_inland() does not exist".  I want my function to be this:

UPDATE inland SET lt_dispatch_level = obs.named_lt_dispatch_level
FROM obs
WHERE obs.created_at = (SELECT MAX(created_at) FROM obs)
AND inland.cartodb_id = 1

How can I create this function and subsequent trigger in CartoDB?  I have been reading the documentation for both Postgres triggers and functions but I need some help.


texts:
Trying to manually create a worldfile for gif images of surface currents supplied by the navy here. I have QGIS and am decent at using it but the world file tool is failing without error on me (even after recompiling GDAL).

So I have both been manually calculating world files and using this website.

Repeatedly, I am coming up with the same numbers, but they can not be correct, at least not in the vertical. The resulting world file places the four corners in the correct spot, but pixels are incorrectly distorted in the vertical, and I believe it is only wrong in the vertical-see image of georeferenced TIF vs world file oriented GIF with Gulf lease areas for reference (those boundaries should not be on the coast). I have gone back to square one and georeferenced the source image several times, checked my calculations, and attempted to track down a QGIS expert to help get the tool working.

Can anyone help me generate the correct numbers for this world file?

What gets me is that I've been able to correctly do this for other images from the NRL site.... this one is just always slightly off.

Current (almost correct but at least line 4 is wrong) world file:
1) 0.02271210592255126
2) 0.00000
3) 0.00000
4) -0.02039280161507402
5) -98.38164194703873
6) 33.08804759959623

GCP file from QGIS if anyone wants a look: https://drive.google.com/file/d/0Bzky2_FkVtI9LXdrbjlENkc2N3M/view?usp=sharing


texts: Existing web tools to digitize non-georeferenced map images?

texts:
Are there any open source web based solutions that allow you to upload an image, georeference it to a baselayer, then draw and export geojson feature data?

Possibly similar to http://geojson.io  .
But includes the ability to upload images and match them up to a baselayer.  And then digitize the features with a drawing tool.

Seems would be possible to create with Leaflet plugins like   Leaflet.DistortableImage

Wondering if someone has already created something similar?


texts: Stretch to MinMax results in blank render output in QGIS?

texts:
I am a new QGIS user, and I can't figure out if I encountered a bug.
If I have a multi-layer image, and "Multiband color" -> "Stretch to MinMax" selected I get the following behaviour. QGIS takes a bit of time to render and the output is blank. I need to force refresh of the view-port, and out of 5 renders randomly 1 time it renders correctly.

If I select "Stretch and clip to MinMax" then everything is fine, except I get clipping.

Blank View-port looks as follows:


texts: Calculation for conversion of Lat/Lon to World Geographic Reference System (GEOREF)

texts: Does any one have a reference to, or the actual code, to convert Lat / Lons to Georef code (http://en.wikipedia.org/wiki/Georef)?  I am also interested in finding a conversion from lat/lon to MGRS.


texts: What does ArcGIS Server error (The required port '4000 ' is already in use) indicate?

texts:
I tried to publish my map in ArcGIS for Server 10.1, but when doing so the following error appeared:

The required port '4000 'is already in use.'

For your information I use ArcGIS for Desktop 10.1.


texts: Seeing Attachments from File Geodatabase at ArcGIS Online?

texts:
I have a feature class (Buildings) related to a table in a file geodatabase (Assets) (1-M RC). The Assets table has attachments of its own (PDFs): mean it has another table in RC with blob data that contain the PDFs.

The Building feature class is in a feature dataset and both relationship classes are in the feature dataset.

I have tried to publish the GDB a few ways (As A full GDB and as a feature service) into AGOL but I can't see the PDFs Attachments. Though I can see my related assets rows while pop-uping the feature's layer entities, I can not  further related it to the attachment table and its contents.

I made sure that the table has its attachments enabled through the REST platform.

My question is this:
1. How to see my PDF's using pop-ups interface in AGOL when interacting with the feature layer (Buildings)?
2. Alternatively, how to open attachments that is related to a stand-alone table in AGOL?

Further question:
1. Where does the attachments table is saved in the AGOL server? Why can't I see it?


texts: Python script for identifying duplicate records (follow up)

texts:
I am trying to calculate a field that marks duplicates and sequentially numbers them.  I am using the code below that I found here: Generating sequential numbering for duplicate values in a field sorted by ObjectID.

It works, but I get a duplicate starting value before it begins sequentially numbering unless I only run the field calculator on the selected values.  When I calculate after manually selecting the duplicate rows and running it on only those values, it works fine.

Here's the code:

prevFieldValue = ''
counter = 0
def GetDuplicateCounter(myFieldValue):
  global prevFieldValue
  global counter
  if myFieldValue == prevFieldValue:
    counter += 1
  else:
    counter = 1
    prevFieldValue = myFieldValue
  return counter

Results running it on all roads (nothing selected):

Results running it after manually selecting the duplicated fields (not too efficient!).


texts: Affine transformation of different sized polygons layer to equal area units

texts: I would like to transform a layer of polygons of different sizearea to become equal in areasize. I have QGIS 2.18 and GRASS GIS 7 available and some basic Python knowledge. Does anyone have an idea of how to do that?


texts: How to use getfeatureinfo with MapServer and OpenLayers?

texts:
I'm trying to use the getfeatureinfo function of Openlayers with MapServer. When I click on a point on my map, I get the following error in Firebug:

<?xml version='1.0' encoding="ISO-8859-1" standalone="no" ?>
<!DOCTYPE ServiceExceptionReport SYSTEM "http://schemas.opengis.net/wms/1.1.1/exception_1_1_1.dtd">
<ServiceExceptionReport version="1.1.1">
<ServiceException code="LayerNotQueryable">
msWMSFeatureInfo(): WMS server error. Requested layer(s) are not queryable.
</ServiceException>
</ServiceExceptionReport>

I've included the javascript and mapfile that I am using. As far as I can tell, I have the correct arguments in my Mapfile but it does not appear to enable querying using getfeatureinfo. I'm also a little unclear on how my html templates should be structured or if they are essential for this.

This is the javascript (I am using OpenLayers 2.12):

    info = new OpenLayers.Control.WMSGetFeatureInfo({
    url: 'http://mysite.com/mapserv.cgi?MAP=points.map',
    title: 'Identify features by clicking',
    queryVisible: true,
    eventListeners: {
        getfeatureinfo: function(event) {
           console.log(event.text)
        }
    }
});
map.addControl(info);
info.activate();

And this is my mapfile, points.map:

MAP
  IMAGETYPE      PNG
  EXTENT         -9780000 2770000 -8850000 3650000
  SIZE           500 600
  SHAPEPATH      "/home/data"
  IMAGECOLOR     255 255 255
  PROJECTION
    "init=epsg:900913"
  END
  CONFIG          "MS_ERRORFILE" "/home/tmp/ms_error.txt"
  DEBUG 5
  WEB
    METADATA
        "wms_title"           "WMS Demo Server"
        "wms_onlineresource"  "http://mysite.com/mapserv.cgi?MAP=points.map&"
        "wms_srs"             "EPSG:900913"
        "wms_enable_request"  "*"
        "ows_enable_request"  "*"
    END
  END

  ##### symbol definition
  SYMBOL
    NAME "circle"
    TYPE ellipse
    FILLED true
    POINTS
        1 1
    END
  END

  ####### layer
  LAYER
    NAME            points
    CONNECTIONTYPE  POSTGIS
    CONNECTION      "host=localhost dbname=data user=postgres port=5432"
    DATA            "geom from points using unique gid using SRID=3086"
    HEADER          'head.html'
    FOOTER          'foot.html'
    METADATA
        "wms_enable_request" "*"
        "gml_include_items" "all"
        "wms_include_items" "all"
    END
    STATUS          ON
    TYPE            POINT
    PROJECTION
        "init=epsg:3086"
    END
    CLASS
      NAME 'point symbols'
      STYLE
        SYMBOL "circle"
        COLOR         100 100 100
        OUTLINECOLOR  10  10  10
        SIZE 50
      END
    END
  END
END

My template files are very basic - just <html><body> for the HEADER and </body></html> for the FOOTER.

Where am I going wrong in this overall process?


texts: Doesn’t cal AGSSOM program to manage AGS mapservices from python script

texts:
We are using ArcGISServer 10.0 on a windows server. Because of database management we need to stop several ArcGISServer Services. For this purpose I use the AGSSOM program, which is called from the following Python program:

import os, subprocess

def manageMapServices(command, svc, svcType):
    PathToAGSSOM = r"c:\Program Files\ArcGIS\Server10.0\AGSSOM_GESTAO\AGSSOM"
    server = "AGSSOM"
    p = subprocess.Popen(PathToAGSSOM + " " + server + " " + command + " " + svc + " " + svcType)

# Calling AGSSOM - stop servico
manageMapServices("-x", "GeoPlaneamento/AMAE2001", "MapServer")

When calling this program I receive the following error:

*Could not find an available ArcGIS Engine or Desktop product to initialize. Could not connect to server. Make sure this user is a member of the AGSADMIN group. Details: The RPC server is unavailable. (Exception from HRESULT: 0x800706BA) AGSSOM [usage]: {server} [-s|start, -x|stop, -r|restart -list|list_services
-listtypes|list_server_types, -describe|describe],-publish|publish from MX D path as second parameter [servicename (or all)] {servicetype|default:MapServer, GeocodeServer, GPServer, GlobeServer}*

This message points to problems with permissions, however the user has all the correct permissions. And when trying to call the AGSSOM program directly from the command prompt I do not receive any error.

Does this anyone know what might be the problem causing this error?


texts: How to add table of non-geo data to Web AppBuilder for ArcGIS?

texts:
I have a parcels layer for a city, and every parcel has ID and I have a table in the database (not linked) contains transactions that happened on these parcels like (buying or renting) and some other data.
Every transaction has the ID of that parcel.

My question is how can I upload this table (transactions) to my Arcgis online map and link every transaction to the parcel using the parcel ID .


texts: Text xyz coordinate import gives data that is not usable in QGIS

texts:
I am using QGIS to try and create a raster from 5m XYZ text point data. I imported using the tab delimited data function and points were created and located where I want them to be. I now want to convert them into raster but it appears the program is not recognizing my Z values. I could not find the Set Z Value function anywhere under Processing and tried the Vector to 3D function under the GRASS Processing menu. However it said the data source was not valid and I got no output. I wonder what might I have done wrong?
I saved the point file as shapefile but it still appears with a strage icon, the three dot icon attached.

Here is a snip of the data which looks ok in the table and functions as a point layer until I try to do any height-related calculations.


texts: Using a custom XYZ tile server in Manifold

texts:
I've a custom tile server that returns a tile based on the X/Y/Z values provided in the URL. The thing is just brilliant when used as a gMaps layer.

I can simply add a new Tile Server to QGIS and call it a day, but now a client of ours would like to import the layer to Manifold which doesn't have support for that (I think?).

I've tried exporting the tiles from QGIS using GeoTIFF but the process is way too slow and I failed to find a way of importing the resulting vrt file to Manifold, plus the resulting images had a ton of weird artifacts.

Ideally I'd share the original data with them, but that is not an option.

Should I use a different format when exporting?

Stitch the ~2m tiles and just have a single image spanning the entire US?

Would Radian be of any help there?


texts: GEE - S1 imagery appears as RGB instead of grayscale

texts:
I would like to visualize a collection of VV Sentinel 1 imagery in Google Earth Engine. This is what I have done so far:

<// Filter the collection - VV only depending on geometry and date
var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
    .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
    .filterBounds(geometry)
    .filterDate('2019-10-10','2019-10-14')
    .select(['VV']);
print(s1);

//visualisation parameters

 var viz = {
  band: 'VV',
  min: -50.0,
  max: -10,
};

// display each image in collection
function addImage(s1) {
  var id = s1.id;
  var scene = ee.Image(s1.id);
  Map.addLayer(scene, viz, id);
  Map.centerObject(geometry,9);
}

// use map on client-side
s1.evaluate(function(s1) {
  s1.features.map(addImage);
});>

Why do the results visualize as RGB? And how to avoid that?


texts: Converting raster layer to polygon form using ArcMap?

texts:
I am trying to convert a raster layer of U.S. climatic regions into polygon format using ArcMap 9.3 under Windows XP. The original data are to be found at the following address:

Koppen Climate Classification from Idaho's Geospatial Data Clearinghouse

Initial attempts to use the Raster to Polygon tool generated "999998 Unexpected Error" and even "Serious Application Error" on a couple of occasions.

To investigate further, I then replicated the procedure on another machine (running version 10 under a more recent operating system), which threw up the more meaningful error: "Raster must be integer type".

Having tried the "Math - Int" tool, it elicited a "000367 Invalid GP data type” error, so I resorted to  the "Math - Times"  tool, successively trying multiples from 100 to 1000000, every time receiving the error: "000864 The input is not within the defined domain”.

Can anyone suggest a workable procedure that I might use to convert the raster layer of the source data into a polygon format?


texts: Find and replace text based on feature class attribute

texts:
I've been struggling with this problem for a few days. It seems to be more programmatic than geospatial in nature, so I'm open to suggestions to redirect this elsewhere.

I have a text file representing a project file for another program. I'm trying to search that text file based on a feature class (or shapefile) with an ID field (ATTR_1) with elevation information (ElevMeters) and replace default terrain heights with heights pulled from the feature. I can't really search for the ID info alone because it appears numerous times throughout the project file. From my perspective, I have to search through the project to find lines with $TERRAIN_HEIGHT :, look three lines above to match $PHOTO_NUM with the equivalent ATTR_1, and then replace the $TERRAIN_HEIGHT value with ElevMeters.

I've gotten my script to property assign and print correct ElevMeters values for each ATTR_1 found. I've been unable to get it to replace the $TERRAIN_HEIGHT : line properly. It prints all new values at the end of the document rather than the line I thought I was directing it to.

import os, arcpy

pdir = r"input directory path"
infile = r"input project file path"
infc = r"input feature class"
fields = ['ATTR_1', 'ElevMeters']

fpath = os.path.join(pdir, infile)
f = open(fpath, 'r+')
searchlines = f.readlines()

for i, line in enumerate(searchlines):
    if "$TERRAIN_HEIGHT : " in line:
        print i #prints line number
        for l in searchlines[i-3:i-2:2]: #go up three lines, skip every other character
            print l #prints photo_num
            print line #prints original terrain_height
            with arcpy.da.SearchCursor(infc, fields) as cursor:
                for row in cursor:
                    if str(row[0]) in l:
                        newline = str(line)[:20] + str(row[1])
                        f.write(line.replace(line, newline))
f.close()

To prove that I'm grabbing values correctly, when I run this code:

for i, line in enumerate(searchlines):
if "$TERRAIN_HEIGHT : " in line:
    print i,
    for l in searchlines[i-3:i-2:2]:
        print l,
        print line,
        with arcpy.da.SearchCursor(infc, fields) as cursor:
            for row in cursor:
                if str(row[0]) in l:
                    newline = str(line)[:20] + str(row[1])
        print newline

I get an output that looks something like this, with line number, photo number, a default $TERRAIN_HEIGHT : of 650.000000 and correct new values underneath:

...
4153   $PHOTO_NUM : 60858
  $TERRAIN_HEIGHT : 650.000000
  $TERRAIN_HEIGHT : 446.301
4165   $PHOTO_NUM : 60859
  $TERRAIN_HEIGHT : 650.000000
  $TERRAIN_HEIGHT : 447.882
4181   $PHOTO_NUM : 60860
  $TERRAIN_HEIGHT : 650.000000
  $TERRAIN_HEIGHT : 448.934
4197   $PHOTO_NUM : 60861
  $TERRAIN_HEIGHT : 650.000000
  $TERRAIN_HEIGHT : 441.156
...


texts: Making shapefile by tracing over Google Earth image?

texts:
I'm trying to create a very simple polygon shapefile of some new areas defined by boundaries specified by my company.

These boundaries don't match any known shapefile as it's a mapping of their infrastructure. They've given me a hardcopy scan of the boundaries traced over a google map.

What I need to do is create a simple shapefile from this, bring into R and match coordinates of another dataset to each specific polygon so I can run some analysis in R based per region.

What would be the quickest and easiest way to do this?


texts: extract() in R raster package slow with non-native format

texts:
I've noticed that running extract takes very different amounts of time depending on whether or not the raster is in memory, on disk in the Native format (.grd), or on disk in another format (e.g., GeoTiff). The following MWE demonstrates:

library(rgdal)
library(raster)
# Make a big raster
rst <- raster(ncol=10000, nrow=10000, xmn=-1000, xmx=1000, ymn=-100, ymx=900)
values(rst) <- runif(ncell(rst))

# Make a big set of points
ext <- extent(rst)
points <- data.frame(lng=runif(100000, ext@xmin, ext@xmax),
                     lat=runif(100000, ext@ymin, ext@ymax))

# Save raster to disk as grd (Native) and tif (GeoTIff)
writeRaster(rst, "tmp_rast.grd")
writeRaster(rst, "tmp_rast.tif")

# Load rasters from disc
rst.grd <- raster("tmp_rast.grd")
rst.tif <- raster("tmp_rast.tif")

# Test extract time for rasters in memory, from native format, and from tif
system.time(test <- extract(rst, points)) # <1s on my machine
system.time(test <- extract(rst.grd, points)) # 6s on my machine
system.time(test <- extract(rst.tif, points)) # 88s on my machine

This isn't too much of a surprise, given that the raster docs note that the Native format is binary. But I'm wondering if there are any useful workarounds here if the input file is in GeoTIFF format (outputted from gdalwarp), assuming the raster is too big to load in memory. Is there a fast way to convert it to binary on the fly?


texts: How to select different values from two columns through Select By Attributes?

texts:
I'm using ArcMap 10.3. I'm trying to select by attributes different values from two columns.

I want to define the areas where

Main rock is GRAVEL And
Sub-rocks are :

a. mudstone tephra loess

b- sand mud

c- ...

d- ...

e- sand mud peat

I have tried to do it as shown below

"MAIN_ROCK" = 'gravel'  AND "SUB_ROCKS" = 'sandstone siltstone mudstone limestone silt lignite' AND "SUB_ROCKS" = 'mudstone tephra loess'  AND "SUB_ROCKS"= 'sandstone mudstone tephra silt lignite' AND "SUB_ROCKS" = 'sand mud peat' AND "SUB_ROCKS" ='sand mud'

After I clicked on apply, nothing was highlighted.

Any idea how can I do this?


texts: Extracting Vertices / Snapping to points from Vector Tile Service (VTS) Layer

texts:
I am using a new vector tile service for background mapping.
I can select/interrogate individual objects from the vector layer with the info tool - but they're not recognised by the snapping tool.
I wanted to be able to snap to the vertices from the VTS layer when creating my own layers on top.
Any ideas/suggestions very welcome!


texts: ArcMap - wrong/no features selected when using select tools

texts:
All of a sudden my features aren't selecting properly when using the "select by rectangle" tool. Previously, I could click on an individual point feature and that feature would be selected.  If I enclosed a group of points with the rectangle, each of these point features would be selected (as expected).  For some reason, this selection behavior is not working.  It is as if each point feature has an invisible line extending some distance to the right, and if that line intersects the "select by" rectangle, then the point feature is included in the selection.

This behavior is consistent for all feature classes, regardless of geometry.

If it helps. . .
Many of my feature classes have heavily edited tables (field edits) and geometries (point features moved)
I am using ArcGIS for Desktop v10.6 (Windows 10)
My data is all stored in a file geodatabase.


texts:
I'm trying to run this very simple script to check a problem in another more comlpex one.

import arcpy
arcpy.env.workspace= 'C:\\Users\\Umberto\\Desktop\\Prove_Script\\Karakoram\\temporary_01'

rst= arcpy.listRasters ()
print rst

I didn't find out what's the matter so far. In python console within arcgis it works without any problem. It gives me this error:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.1\Lib\site-  packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__
  File     "C:\Users\Umberto\Documents\UniversitᝐhD\Karakoram\MODIS_Pakistan\Script\test_listClip.py",     line 4, in <module>
    rst= arcpy.listRasters ()
AttributeError: 'module' object has no attribute 'listRasters'


texts: Correct area calculations: When to use planimetric, when ellipsoidal?

texts:
Sometimes I am not sure whether to use planimetric or ellipsoidal model to calculate areas of polygons. If there is a "reference" data (e.g. a given shape_area column) I compare results and take the given one.
But is rule (of thumb) when to use planimetric and when a ellipsoidal model?


texts:
I need to generate hundreds of MBTiles from different shapefiles and same styles, so I created a xml from tilemill with a single shapefile and styling, using python landez library to generate MBTiles from this xml (keep changing datasource for each shapefile). Problem is landez is taking too long to create while tilemill is faster in exporting MBTiles, also through landez, few tiles are missing at different zoom levels (blank tile, so default color).

I was just wondering if this can be done using tilemill, (tilemill seems to export only projects).

PS: I want to do this process through terminal.


texts:
We're making a light web front to a database, and I saw this in a test screenshot:

Despite those coordinates being nowhere near North America, the "latitude/longitude" label for those coordinates really ruffled my jimmies. This got me wondering if there is a generic yet unambiguous term for each axis.

The first idea that came up was to just change the labels to X and Y Easting and Northing, but my understanding is those are terms used specifically for Cartesian coordinates. If the CRS was 4269, then the labels we have now would be correct.

We can probably derive the proper axes labels from the CRS, but I thought I'd ask if there's an accepted generic term (other than "value that increments along the north axis") to just sidestep the issue.


texts: Shapely .project() function not working on horizontal or vertical sections of line

texts:
I have a function called intersect_point_to_line that includes the following relevant code snippet:

for i in range(len(segment.coords[:])):
    if i < len(segment.coords)-1:
        seg = LineString([Point(segment.coords[i]),Point(segment.coords[i+1])])
        if seg.distance(point) <= d:
            intersect = seg.interpolate(seg.project(point))

This is designed to take a list of points that are slightly off of the line and snap a location that is used to write a new vertex that is on the line.

I have a very large list of points that I am running against a complex LineString that twists around quite a bit. The function works correctly except where the line goes vertical or horizontal. In those sections there is a short area where no points are snapped.

I am confident that the remainder of the code is correct. Am I missing something about the .interpolate() or .project() methods?

Edit: Here is the whole function

def intersect_point_to_line(segment, point, d):
    # iterate through each adjacent pair of points in the line segment
    for i in range(len(segment.coords[:])):
        if i < len(segment.coords)-1:
           seg=LineString([Point(segment.coords[i]),Point(segment.coords[i+1])])
            if seg.distance(point) <= d:
                intersect = seg.interpolate(seg.project(point))
                if not intersect.is_empty:
                    # calculate corrected distance to the actual point of intersection
                    real_dist = point.distance(intersect)
                    intX = intersect.coords[0][0]
                    intY = intersect.coords[0][1]
                    return (real_dist, intX, intY)
                else:
                    return False

The arguments are:
segment - a Shapely LineString of lat/lon pairs in Mercator 3785
point - a Shapely Point in Mercator 3785
d - a decimal value representing a threshold distance

This works for all points tested, except for where the LineString goes horizontal or vertical, for a distance of roughly 30 - 60 feet. I am attaching a representative screenshot below.


texts: Ways to access and download the metadata of OpenStreetMaps

texts: Does anyone know if there are any good methods or apps to access and download the metadata of OSM from desktop?


texts: Is there an equivalent to the Dissolve all function in GDAL dissolve user interface?

texts:
I am looking to compare the speed of the dissolve tool because the QGIS dissolve tool is performing very slow for me. However, I need to dissolve all the features into one so that I can create an outline of the perimeter of these.

The QGIS dissolve tool has this as a tick option (please see first image below).

Is it possible to do this with the GDAL UI?


texts: Serve up an ESRI Tile Cache without a service?

texts:
Trying to take an exploded ESRI tile cache and serve the PNG images on a cloud storage space.  I want to be able to generate the tile caches from ArcMap (from rasters and mosaic datasets), but serve them as tiles that don't require an ArcGIS Server service. Problem is there will be hundreds of these rasters each requiring a separate service.
It seems that the folder structure that ESRI uses is different than the standard slippy map structure. Not sure if there is any documentation on it.
http://[abc].tile.openstreetmap.org/zoom/x/y.png
https://{subDomain}.tile.openstreetmap.org/{level}/{row}/{col}.png

Example of the ESRI tile cache structure. The row folders aren't named with simple integers like in OpenStreetMap. Seems to use hexadecimal numbers instead.
L10
    R0000015c
        C000000b0.png
        C000000b1.png
        C000000b2.png
    R0000015d
L11
L12
L13
L14
L15
L16

The WebTileLayer from the ArcGIS API for JavaScript doesn't work with non-ArcGIS Server map files.
https://developers.arcgis.com/javascript/latest/api-reference/esri-layers-WebTileLayer.html
How can tiles from an ESRI TileCache be served up just as a folder of images without a service?


texts: Pyqgis Set the field value of a dynamically created vector layer in QGIS 3 before user saves to database

texts:
I am quite new to both QGIS and Python.
I have inherited a project that needs to be upgraded from using QGIS 2.0 DUFOUR to 3.8 Zanzibar.

The way it works is, there is a separate application written in asp.net. On one of the web pages there is a button called 'Digitize'.

When the user clicks this it saves a file with the details of the current record, onto the web server.
The user then opens QGIS via a shortcut on their desktop. The QGIS app is locked down with only minimal functionality.
When the project has finished loading, a plugin button becomes available on the toolbar. This button creates a new Vector Layer by connecting to a mssql database.

The vector layer name is based on dynamic values in the file on the server, so it changes each time.
The layer is also removed every time QGIS is started.
It displays the features/polygons that have already been added (if there are any).

The user clicks the pencil to begin editing.
The user draws a polygon.
The user clicks save.
The user is always prompted to enter an id which the record needs as a foreign key in the database. However, the user will never have access to this ID which is a long string e.g. '0000000001-0000000000000264100' that is hidden in the code.

I can get this value in the code from the server file (see sublayerref below).
My problem is I cannot find a way to set this programmatically.

Does anybody know how to do this?

Here is the code so far, the best I can come up with at the moment is using the signal 'editingStarted', but I haven't found any success setting the value of 'ForeignKeyFieldId' field.

    {
    def feedback2():
    print("XXXFEEDBACK2 ACTIVATED:" + sublayerref)
    df = QgsDefaultValue('var("varsublayerref")', True)
    for field in vlayer.fields():
        print("LayerField:" + field.name())
        if field.name() == "ForeignKeyFieldId":
            field.setDefaultValueDefinition(df)
            print("Job done")
crs = QgsCoordinateReferenceSystem("EPSG:27700")
uri = QgsDataSourceUri()
uri.setConnection(servername, portnumber, databasename, logon, password)
uri.setDataSource( schemaname, tablename, geometryfieldname,  sublayerid + " = '" + sublayerref + "'" )
uri.setSrid("EPSG: 27700")
uri.setUseEstimatedMetadata(True)
uri.setWkbType(QgsWkbTypes.Polygon)
primaryKeyColumn = "RecordId"
uri.setKeyColumn(primaryKeyColumn)

symbol = QgsFillSymbol.createSimple({'line_style': 'dash', 'color': layercolour })
symbol.setOpacity(0.5)
vlayer = QgsVectorLayer()
vlayer.setCrs(crs)
vlayer.setCustomProperty("labeling", "pal")
vlayer.setCustomProperty("labeling/enabled", "true")
vlayer.setCustomProperty("labeling/isExpression", True)
vlayer.setCustomProperty("labeling/fontFamily", "Arial")
vlayer.setCustomProperty("labeling/fontSize", "10")
vlayer.setCustomProperty("labeling/fontWeight", "2")
vlayer.setCustomProperty("labeling/fontBold", True)
vlayer.setCustomProperty("labeling/placement", "1")
vlayer.setCustomProperty("labeling/fieldName", "'" + labelprefix + " ' || " + '"' + pkfieldname + '"')
vlayer.setOpacity(0.5)
vlayer.updateExtents()
vlayer = iface.addVectorLayer(uri.uri(), layername, 'mssql')
vlayer.renderer().setSymbol(symbol)
vlayer = iface.activeLayer()
vlayer.editingStarted.connect(feedback2) ##testcode
}


texts: Using Click locations to determine azimuth in QGIS 3

texts:
I am using QGIS 3 64bit. While adding points to a shapefile (layer), what expression captures the insertion point coordinates of a point feature to use as a start point for an azimuth determination followed by a second click along a known trend (eg. from a georeferenced layer image) the coordinates of which are used as the end point and use these start and end points to determine (and enter in the shapefile table) the azimuth bewteen them?

My shapefile table has columns id, station, azimuth. The azimuth value is used to orient the point symbol using a rule-based command.


texts: Performing localized raster edits using ArcGIS Spatial Analyst?

texts: I have a raster image of water depth, float point values 1-15, in a very small area some erroneous values are in the raster with values of 100 - 375.  I'm looking for a repeatable work flow in ArcMap 10.5, Spatial Analyst level.  I've used one work flow to set the values to nodata, one work flow to set to a fixed float point value, but what I'd really like is to create polygon or raster mask, and then have the grids of the original raster, within the mask area, to be replaced with nearest neighbor values from the original raster outside the mask area.


texts: QGIS Essen ruler problems with EPSG:3857

texts:
I am using QGIS version 2.14.3-Essen. I am pretty new to this so bear with me. I have 3 layers: one of river polygons in MN, one of lakes, and the last an OpenLayers plugin of Google Satellite.

I am using the ruler, and originally 90% of the data matched with previous information we had (I am double checking, you see). However, I lost my project and had to re-install the 3 layers. Now the ruler consistently measures about 1.5 times what it should be. For example, the distance should be 594 meters away, but when I use the QGIS ruler it now says it is 892 meters.

I thought it could have been the order in which I download the layers, but that did not seem to solve it. Then I thought a prj issue, however I cannot add the prj files separately as vectors. Other sources says that the prj is not the issue because the shp file will automatically read the prj (as well as the other files). So then, I thought it could be a "Set Layers CRS" problem. I am now in EPSG: 3857, but when I find the CRS for my area (Minnesota EPSG: 26793) it layers my 2 water body files over Africa...

This person seems to have had a similar problem:
http://lists.osgeo.org/pipermail/qgis-user/2010-December/010597.html

Their solution was to: "get correct measures by reprojecting your shapefiles into a suitable UTM srs [and] may have to change the project crs too." However, I am not exactly sure what this entails..

Any advice would be great!


texts: Not able to add markers on the OpenLayers map

texts:
I have a following short HTML code, which displays a map. I can't find any way to add a marker on this map. I've tried different examples, but everything's failed so far. Here I post one of the codes I've tried.

    <link rel="stylesheet" href="https://cdn.rawgit.com/openlayers/openlayers.github.io/master/en/v5.1.3/css/ol.css">


texts: CenterPoint of Polygon in ArcEngine .Net 10

texts:
I am currently trying to get the center point of a polygon. In GIS, i can see the "x" mark in the center of the polygon, but I do not know how to get it from the IPolygon object.

I tried getting the enevelope and getting the midpoint from that but that didn't work out very well, as getting the envelope does not always ensure to get the midpoitn where marked "x".

Thanks and Regards,
Kev


texts: Problems viewing ASTER imagery in google earth engine

texts:
I am trying to view ASTER imagery for Mount Spurr from 2004 to 2006.
I am loading the imagery in a list "listOfImages" and access them for viewing via listOfImages.get(). Despite all images appear in my list and are loaded as layers, only a fraction of them are visible on the map when clicking on the specific layer. The code I use is as follows:

// Load a landsat image and select three bands for Mount Spurr
// Mount Spurr
// start: 2004-03-01
// end: 2006-05-30

 var dataset = ee.ImageCollection('ASTER/AST_L1T_003')
              .filterDate('2004-03-01','2006-05-30')
              .sort('system:time_start',false)

// filter imagery
var Image_filtered = dataset.filterBounds(geometry_Mt_Spurr)
 print(Image_filtered.size())

// convert to list
var listOfImages = Image_filtered.toList(Image_filtered.size());
print('List:',listOfImages);

// Get the size of the image list (this is a server side ee.Number object).
var len = listOfImages.size();
// build and show layers in GEE
len.evaluate(function(l) {
  for (var i=0; i < l; i++) {
    var img = ee.Image(listOfImages.get(i));
    Map.addLayer(img, {bands:['B3N','B02','B01']}, 'Img'+i.toString(),false);
    print('Img-'+i.toString());
    print(img.date());
  }
})

The ASTER imagery has been loaded as an ImageCollection (ee.ImageCollection) and geometry_Mt_Spurr is just a rectangle over the area of interest.

Working with less images in the list did not solve the problem.


texts: Open Source Remote Sensing vs eCognition, ERDAS etc Is it worth it?

texts:
What, if any, would be the main benefit of using a commercial tool like eCognition (in particular) or ERDAS as opposed to something like Orfeo Toolbox or OpenCV? (Feel free to suggest other open source software for RS if you think it fills a gap also.)

Put aside the GUI as I would likely be scripting it anyway.  It seems to me that most of the fundamentals of Remote Sensing are available as easily scriptable OSS elsewhere.

Things like mosaicing, edge detection, basic band calcs e.g. (Band_1 - Band_2) that can be strung together for things like NDVI, supervised and unsupervised classification, high/low pass filters etc all seem to be available.  Furthermore the basic C++ programs written for these are no more complicated than Python given the level at which the methods are exposed to the user.

I am assuming the commercialproducts shine when it comes to more advanced stuff like object recognition and the like.  Also I assume they would have the benefit of a helpdesk or support if required.

I have used the free trial for eCognition but given that I am usually doing pretty basic stuff I never get a chance to see where it may be better.  I am also assuming that it takes a while cutting your teeth on the basics in a program to realise the helpful time savers and get the most benefit out it.  Time you don't really get when using a trial version.

Is there anything I am missing?  Anything that completely stands out that you are thinking 'Wow! this guy is an idiot for not knowing about capability X'? Conversely do you think it is true that most of this stuff is basic technology that has been around for years so it is only the more advanced part of the market pushing the technological envlope that requires these tools?


texts: Mapbox expression to return greater of two values

texts:
I'm looking to plot the equivalent of the SQL query

SELECT column_3
FROM my_object
WHERE map_layer is not null
AND column_1 > column_2;

Using this answer I can see where the expression should go (in the filter) and that I can combine multiple filters but not how to compare values.

Which gives me

"source-layer": "my_object",
        "filter": ['all',["==", ["typeof", ["get", "col_3"]], "string"],[another expression here]]


texts:
I have an excel table. I created a database table in SQL Server database platform in ArcMap. I added fields, fields' data type, and domains to some data fields (columns). When I open the table, all I see is just a table with field names only.

Can I not view the database table in SQL Server in ArcGIS desktop?

Also, how can I populate rest of my records ( rows) in the table?