texts:
Is it possible to scale each polygon by a numeric value in an attribute table field using arcpy. I have updated the code from some previous excellent work done by Evil Genius

import arcpy
import math

def scale_geom(geom, scale, reference=None):
   """Returns geom scaled to scale %"""
   if geom is None: return None
   if reference is None:
    # we'll use the centroid if no reference point is given
    reference = geom.centroid

refgeom = arcpy.PointGeometry(reference)
newparts = []
for pind in range(geom.partCount):
    part = geom.getPart(pind)
    newpart = []
    for ptind in range(part.count):
        apnt = part.getObject(ptind)
        if apnt is None:
            # polygon boundaries and holes are all returned in the same part.
            # A null point separates each ring, so just pass it on to
            # preserve the holes.
            newpart.append(apnt)
            continue
        bdist = refgeom.distanceTo(apnt)

        bpnt = arcpy.Point(reference.X + bdist, reference.Y)
        adist = refgeom.distanceTo(bpnt)
        cdist = arcpy.PointGeometry(apnt).distanceTo(bpnt)

        # Law of Cosines, angle of C given lengths of a, b and c
        angle = math.acos((adist**2 + bdist**2 - cdist**2) / (2 * adist * bdist))

        scaledist = bdist * scale

        # If the point is below the reference point then our angle
        # is actually negative
        if apnt.Y < reference.Y: angle = angle * -1

        # Create a new point that is scaledist from the origin
        # along the x axis. Rotate that point the same amount
        # as the original then translate it to the reference point
        scalex = scaledist * math.cos(angle) + reference.X
        scaley = scaledist * math.sin(angle) + reference.Y

        newpart.append(arcpy.Point(scalex, scaley))
    newparts.append(newpart)

return arcpy.Geometry(geom.type, arcpy.Array(newparts), geom.spatialReference)

incur = arcpy.da.SearchCursor(r'D:\ZOC\POLYPGDIS_MASTER.shp', ['OID@','SHAPE@'])
outcur = arcpy.da.InsertCursor(r'D:\ZOC\POLYPGDIS_MASTER_1.shp', ['SHAPE@'])

for row in incur:
    # Scale each feature by 0.5 and insert into dest_fc
    outcur.insertRow([scale_geom(row[1], arcpy.da.SearchCursor(r'D:\ZOC\POLYPGDIS_MASTER.shp', ['Scale_Fact']))])
del incur
del outcur


texts:
I am using TileCache from MetaCarta to create tiles for a WMS service. It creates a tile with an URL like
http://example.com/tiles/Admin/06/000/000/357/000/000/158.png.

I am trying to understand how to start from Map Coordinates, and get the appropriate tile.

I am looking for some documentation which explains the Folder structure/hierarchy of the server cache directory


texts: Calculating polygon main direction and creating line out of it

texts: I would like to find out the orientation of a polygon and create a line of that direction.See image below.


texts: Issue with speed of adding new features to a layer

texts:
I have a list of features that I would like to add to a layer. The features are built from a list of objects that contain a number of attributes as well as a set of coordinates.

var lines = data; //Data represents the object list

for( var count = 0; count<lines.length; count++)
{

   if (lines[count] != null) {
      var latitude = lines[count].latitude;
      var longitude = lines[count].longitude;

      var pt = new Point(xy[0], xy[1], new SpatialReference({ wkid: app.map.spatialReference.wkid }));
      pt.attributes = {};
      //Code to set each attribute is here
      pt.attributes.attributeName = attributeValue //Taken from the data list for each value - left out all data transfers to preserve space

     var graphic = new esri.Graphic(pt, null, pt.attributes);
     Layer.applyEdits([graphic], null, null);
                            }
                        }
                    }
                }

The process works and each feature is drawn on the map, however the process is super slow.
Accessing the console window throws no errors and ArcGIS Server manager also throws no errors. However checking the network tab in the developer window shows a list of apply edit functions pending, waiting them all out for a list of around 100 can take around an hour which is too long.

Is there some kind of process that is slowing this down? Since debugging my code reveals no issues, the issue seems to be once the apply edit functionality is called for all the features


texts: Shifting extent of overview map, keeping scale, centering on extent indicator in ArcGIS Pro

texts: In ArcGIS Pro 2.7 I've created a map series with a main map frame (1:5000) and an overview map frame (1:250000) with an extent indicator of the main map frame. Once the map series is enabled the extent indicator moves out of my overview map. So the extent indicator should stay in the middle of the overview map frame, the extent of the overview map frame moves. Its like to say "stay focused on the extent indicator, keep scale 1:250 000". If I put the whole map in its full extent in the overview map you wouldn't be able to see the extend indicator of the main map.


texts: How to delete duplicate features from shapefile/OSM XML using code/command line tools only? (NOT QGIS/ArcGis)

texts:
I have shapefile with lines, some of them are perfectly overlaps - duplicate, although each one have unique id etc.

When I view the file in QGIS I see they overlap.
When I convert the file with ogr2osm to osm XML & import to JOSM, and running validation I also get warning of duplicate ways (the lines converted to ways...).

How can I "programmatically" detect and remove such situations? I don't mind doing it in the shapefile or later in the OSM XML.

Note 1 - The attributes (tags) are unique - the actual spatial geometry is duplicate.

Note 2 - I want to use something like ogrinfo, ogr2ogr etc GDAL libraries either as commandline tools or as import library for python/java code. If other code tools exits it's perfect also.


texts: syntax for selectbyexpression MTFCC in pyqgis

texts:
Does anyone know the correct syntax for the selectbyexpression processing algorithm for selecting MTFCC's? This is what I have right now, but it doesn't work so I'm guessing my problem lies in the syntax for the "expression" variable.

layer = qgis.utils.iface.activeLayer()
expression = '"MTFCC" = "H3010"'
processing.runalg("qgis:selectbyexpression", layer, expression, 0)

When I select by expression using the GUI, it shows something like "MTFCC" + "H3010" and gives me an output preview of 'S1400H3010' and I don't understand. Anything else I try just gives me an "expression invalid" error.

I am using pyQGIS v2.14.0


texts: pgr_withpoints() How to select edges geometry till points?

texts:
I have two tables, table poi consist of points with point id pid, nearby/associated edge_id,fraction and point geometry the_geom. And edge table ways which has all the fields required(source, target, cost, reverse_cost, the_geom) for pgRouting.

While trying to get route from poi to poi using pgr_withpoints(), I am bit stuck with getting geometries till the poi. With the query I will get the exact cost till the poi, but I am not able to get the geometries?

select * from pgr_withpoints
('select id, source, target, cost, reverse_cost from ways ORDER BY id',
'select pid, edge_id, fraction, side from poi', -101, -121);

Now to draw lines I can join the pgr results with ways table.

select seq,path_seq,node,edge,pgr.cost,agg_cost,st_astext(the_geom) as geom
from pgr_withpoints('select id, source, target, cost, reverse_cost
from ways ORDER BY id', 'select pid, edge_id, fraction, side
from poi', -101, -121) pgr join ways w on pgr.edge = w.id ORDER BY pgr.path_seq;

But since ways are linestring geometry, route draws till end of the edge. How can I select geometries till poi with pgr_withpoints() ?


texts: Meaning of Simplify's tolerance parameter

texts: The standard parameters of Douglas-Peucker's simplify algorithm are geometry and tolerance (e.g ST_Simplify in PostGIS). What's the meaning of the tolerance parameter? I know that the bigger the value, the coarser the geometry will be. But does the number has any unit or it is just arbitrary?


texts:
I am trying to render cities on the map, using Geonames’ dataset (source). I have downloaded the top 10000 cities as a GeoJSON file, and loaded it into QGIS. However, when I zoom out, this shows too many points: they completely fill the map.
I have enabled labels, and there is the option under Layer Properties › Labels › Rendering › Feature Options to limit the number of features rendered. This way, only the largest cities in an area will be labelled (I have set render priority to be based on population), but when I zoom in I can see more city names appear.
Is it possible to apply a similar limit to the symbols themselves? Can I hide all but the largest 100 cities within the viewport? I do not want to hard-code limits for population and zoom level, because I would like cities in less populous areas to appear sooner when zooming in there.


texts: Parsing field values using ArcGIS Field Calculator?

texts:
I am using 10.0 ArcView and am trying to get the "R" in R12345678910 to drop.

The numbers are in a column in the attribute table and i've always used VB scripting.

I was using Field Calculator.


texts: Resolving ImportError: No module named arcgisscripting?

texts:
I am trying to run an ArcGIS python script from Eclipse and get the following message:

  ImportError: No module named arcgisscripting

I can run this in IDLE... any ideas?


texts: JSON String From VB.NET to Javascript over 100,000 characters

texts:
I have a string of properly formatted JSON from a VB.NET function that I am trying to pass to JavaScript where I can convert it into an OpenLayers vector layer. The string can be variable length, sometimes as small as 5,000 characters and sometimes as long as 40,000,000 characters.

Here is a small example of the JSON:

"{"type":"FeatureCollection","features":[{"type":"Feature","geometry":{"type":"MultiPolygon","coordinates":[[[[-79.546437887619277,41.450810394091974],[-79.546373336415854,41.450878338383419],[-79.546186240728744,41.450930517576978],[-79.546030161587709,41.451004424821683],[-79.546437887619277,41.450810394091974]]]]},"properties":{"tract":11001180000000,"stand_num":"11001180000006","acres":0.585540920000}},.....etc.

If the string is less than 100,000 characters I can successfully send it to my JavaScript function as a normal string variable OR as a serialized variable, can successfully read it as GeoJSON, and add the features to my OpenLayers map using the normal OpenLayers GeoJSON methods. The following JavaScript takes the VB.NET serialized string and adds it to the map;

            var geojson_format = new OpenLayers.Format.GeoJSON();
            var vector_layer = new OpenLayers.Layer.Vector("query_POSTGIS");
            var parsedJSON = JSON.parse(outputStr);
            alert(geojson_format.isValidType(parsedJSON, "FeatureCollection"));

            var geoJsonStr = geojson_format.read(parsedJSON, "FeatureCollection");
            vector_layer.addFeatures(geoJsonStr);

This code works great until I get to the 100,000 character limit. Any suggestions on how to make this work with strings over 100,000 characters?


texts: Swapping bytes for SRTM 16-bit binary format in ArcMap?

texts:
I am using ArcMap 10.2.2.

I have two .hgt files in 16-bit signed integer data. All I have done so far is to mosaic them together (Mosaic to new raster).

Since the DEM data are provided in Motorola or IEEE byte order, which stores the most significant byte first ("big endian"), I`d have to swap the bytes. The only examples I have found show how to swap the unsigned integer data,
e.g. CON([in_grid] >= 32768, [in_grid] - 65536, [in_grid])

Can someone advise me which expression should I use for 16-bit signed integer?


texts: How to make my programmatically added map items show up in the composer items list?

texts:
I try to make a composer via python. I add item, make a composer view and open, but i can't see item in "Items" dock. Also, when i select it, i don't see item properties in "Item properties" dock.

from qgis.core import *
from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *

# make a composition
mapRenderer = iface.mapCanvas().mapRenderer()
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
view = iface.createNewComposer()

# add label
lbl = QgsComposerLabel(c)
lbl.setText("test text")
lbl.setItemPosition(10,10)
lbl.adjustSizeToText()
c.addComposerLabel(lbl)

# view composer
view.setComposition(c)

How i can fix it?

UPD:
I setup model for my item and now i see it in Items dock, but properties still doesn't show....
Additional code:

# find items and items prop widget
win = view.composerWindow()
items_widget = win.findChildren(QDockWidget, 'ItemsDock')[0]
items_prop_widget = win.findChildren(QDockWidget, 'ItemDock')[0]

# set items model
items_list_widget = items_widget.findChildren(QTreeView,'')[0]
items_model = c.itemsModel()
items_list_widget.setModel(items_model)


texts: a localised 2D coordinate system from GPS coordinates

texts:
I spend a few days reading about coordinate systems, map projections, datums, etc. Since I am new to gis, I am not confident I am doing this right.

Problem at hand:

I have a GPS device that gives me GPS coordinates corresponding with the WGS84 datum. That is nice but what I actually want is to have the coordinates in a local (x,y,z) coordinate system with z being the height.

I want to measure distances in an area of about 10x10km. Also I want to pick my own local origin. So I can say the corner of building A at ground level is (0,0,0) and the first floor of building B on the top of the hill is like (5km, 6km, 150m).

I want to be accurate on sub-cm level. The error should be 1cm or less.

Also I want to do the conversion between GPS and local coordinates by myself, in my own code. So I have to understand the underlying mathematics.

Current state:

As far as I understand I would need a conformal projection? Like Mercator or transverse mercator. Web-mercator is apparently a no-no because it's not conformal. A secant projection (instead of tangent) would minimise the deviation in z (height)?

UTM is widely used it seems. But that is for area's much larger than mine, so I maybe I can do better than using a standard UTM zone?

In any case, UTM or a more localised TM projection or something else entirely, I am not sure which workflow and formulas to use.  Some are relatively straightforward, some are rather long like these. In the latter case I also don't see how these were derived/obtained.

Then there is the issue of geodetic and geocentric coordinates. WGS84 (lat, long) is geodetic, but the (x,y,z) are geocentric. Not to mention the issues around height.

The websites I found thus far are either an online converter or a shallow introduction into map projections. Shallow as in little to no math or derivations of formulas.

Can you help?


texts: Multi-select / modify event Openlayers 3

texts:
I'm trying to incorporate a modify feature where the user will move a point and any underlying geometry is also modified. Currently to draw a route on a map I'm taking points from a database and creating a Feature for each(ol.geom.Point) and then pushing the same coordinates into an array for the underlying linestring. All of these are added to a Feature Collection so if there were 28 points in the database, there will be a feature collection with 29 items(28 points and a linestring).

The Select and Modify interactions are set to multi:true and this seems to occasionally allow me to select both the waypoint and the vertex point beneath but it's really hit and miss.

These are sample Modify Start/End functions:

session.interactions.modify.on('modifystart', function (e) {
        session.selectedFeatures = session.activeRoute.Layers[0].getSource().getFeatures();
        session.testArray = e.features.getArray();
        for (var p = 0; p < session.testArray.length; p++) {
             session.modifiedGuid = session.testArray[p].get('Guid');
             session.originalCoordinates = session.testArray[p].getGeometry().getCoordinates();
                }
  });

  session.interactions.modify.on('modifyend', function (e) {
         session.selectedFeatures = session.activeRoute.Layers[0].getSource().getFeatures();
         var lsRoute = session.selectedFeatures[session.selectedFeatures.length - 1];
         console.log(lsRoute);
         console.log(lsRoute.getGeometry().getCoordinates());
         console.log(session.testArray);
                for (var p = 0; p < session.testArray.length; p++) {
                    console.log(session.testArray[p].get('Guid'));
             console.log(session.testArray[p].getGeometry().getCoordinates());
                    session.modifiedCoordinates = session.testArray[p].getGeometry().getCoordinates();
                    updateModifiedGeom(session.modifiedGuid, session.originalCoordinates, session.modifiedCoordinates);
                }
            });

The updateModifiedGeom function was logging to the console the Guid, original, and modified coordinates but since I specified multi: true on the modify interaction it's no longer picking up these details, only the linestring array beneath.

I'm sorry if this seems a little garbled but I've been struggling for months to be able to display rich data points on a linestring and this is the only way I've been able to resolve it to any degree.

I'd previously posted this question which received no responses and I've managed to get this far but hopefully someone can tell me where I'm going wrong. I can't post a fiddle as it would take far too long to adapt the code to work in a stand-alone way.


texts:
I worked with a photo that is apparently a bit offset.
Now I have more then 10,000 points with that offset.
Now that I have the correct map I would like to move them all to their new location.
Is there a way to move them all at once?
I did not manage to do that with the vertex tool on QGIS 3.14 Pi.


texts: How to make sure OnSketchFinished always works?

texts:
In my VBA application I have a listener set up that listens for OnSketchFinished, when the polygon is created a message box appears and then a userform.

My problem is that the listener will work 90% of the time, not 100%.

Here are the snippets of code in ThisDocument:

Private WithEvents EditEvents As Editor
Private m_pEditor As IEditor
Public Sub SetEvents()
    ' This is the listener for the Editor
    Dim pID As New UID
    pID = "esriEditor.Editor"
    Set m_pEditor = Application.FindExtensionByCLSID(pID)
    Set EditEvents = m_pEditor
End Sub

Private Function MxDocument_OpenDocument() As Boolean
    Call SetEvents
End Function

Private Sub EditEvents_OnSketchFinished()
    Call frmCityDev.EditEvents
End Sub

In my form I have code for EditEvents which brings up a msgbox to notify that a new polygon has been created.

Am I missing something to make sure that the listener will always work?


texts: OpenLayers.Layer.XYZ not referencing tiles correctly

texts:
I have tiles stored here:

../CO_054/MtViewBaseMap/MtViewBaseTiles

in the order : C:\Program Files (x86)\GeoServer 2.10.3\webapps\CO_054\MtViewBaseMap\MtViewBaseTiles\0\0\0.png

I'm using this to reference the tiles:

 var tiles = new OpenLayers.Layer.XYZ(
           "CO_054",
           [
            "../CO_054/MtViewBaseMap/MtViewBaseTiles/${z}/${x}/${y}.png"
           ],
           {
            attribution: "Data copyright OpenStreetMap contributors",
            sphericalMercator: true,
            wrapDateLine: true,
            numZoomLevels: 18
           }
         );
        map.addLayer(tiles);

However, I only get broken image icons and the following errors:

  HTTP ERROR 404

Problem accessing /CO_054/MtViewBaseMap/MtViewBaseTiles/11/1025/1021.png. Reason:

    Not Found
Powered by Jetty://

I can see that the directory is not being constructed correctly but it seems that my code is set up to match where the .png files live. I mimicked this from another html file with the same file structure and it works fine.

At CO_054/MtViewBaseMap/MtViewBaseTiles/11/ I have a folder named 329 and two .png files (793.png and 794.png). Why is the WMS looking for 1023/1023.png?

This is what it should be reading:

Here's what the directory looks like in z = 11:


texts: ABC of serving large amounts of data / shapes?

texts:
Some background: web developer, with some experience of Google Maps, currently developing a fairly complicated tool for archeological studies. Outside of Google Maps, I have little to no knowledge of mapping services, gis and whatnot.

I have been looking into OpenData sets to implement, and one which I am particularly interested in raised several questions.

The data is at http://www.datiopen.it/it/opendata/Regione_Lazio_Carta_geologica_scala_1_25000 , tab "Mappa".

My first reaction was, geological data, just the thing, awesome.
Then I saw it runs up to about 16,000 elements, all polys - less awesome, but not yet necessarily fatal.
I finally downloaded the data, and found it is a huge shapefile (35Mb zipped - 55Mb unzipped).
I learned about gdal just about enough to convert the shapefile to an even worse 112Mb kml.

Even stripping away all the extended data from the kml, or picturing it conveniently stored in a database table, it does not look like any amount of data you would wish to serve over an internet connection.
The linked page, however, does it in but a few seconds.
The way it reacts to mouse events makes me think it simply finds out which poly the clicked coordinates fall within, but, here's what puzzles me, at the same time the contours of the clickable areas are already visible on the map before clicking - suggesting the polys are already loaded.

As this is making no sense to me, I suspect Something Everyone Already Knows But I Don't is at work here, possibly starting with the environment.
I have spent some time reading forums and SE itself, but lack of experience leads to lack of perspective.

Can anyone tell by looks what is happening on the linked page, and tell me where to direct my research?


texts: Inverted lat lon on drawn polygons in Openlayers 3

texts:
I am working on a map in openlayers 3 with geoserver I display the layers as WFS and I added the draw and modify interactions, but when I save the drawn polygons the coordinates of the newly drawn polygon are inverted.

Here is my drawing code:

var draw;
var formatwfs = new ol.format.WFS();
function addInteraction() {
    draw = new ol.interaction.Draw({
            features: features,
            type: /** @type {ol.geom.GeometryType} */ ('Polygon')
            });
    // when a new feature has been drawn...
    draw.on('drawend', function(event) {
        var feature = event.feature;
        feature.set('the_geom', event.feature.getGeometry());
        var node = formatwfs.writeTransaction([feature], null, null, {
                        gmlOptions: {srsName: "EPSG:4326"},
                        featureNS: "My_Name_Space",
                        featureType: "My_Name_Space:My_Layer"
                    });
         $.ajax({
                    type: "POST",
                    url: "http://localhost:8080/geoserver/wfs",
                    data: new XMLSerializer().serializeToString(node),
                    contentType: 'text/xml',
                    error: function(e) {
                            var errorMsg = e? (e.status + ' ' + e.statusText) : "";
                            alert('Error saving this feature to GeoServer.

'
                                + errorMsg);
                            },
                        context: this
                    });

                });
        map.addInteraction(draw);
    }

Every Polygon I draw has its vertices coordinates inverted and so it gets into another location.

How can I fix this ?


texts: Show how many Polygons are overlapping in different styles

texts:
First of all I use QGIS. I got a bunch of polygons that describe a 10 min. footway from different addresses.

What I want to do now is to illustrate the areas where the polygons are overlapping. My trouble is how to show the differences between 2 overlapping polygons and 3, 4 or 5 overlapping ones.

My idea is to do some query where I get a Layer with all those overlaps and in one column of the attribute table it says how many polygons are overlapping. I hope the attachment helps to understand what I'm looking for.

Do you have any idea how to do that?


texts: OSM or Google Tile problem fitting to Geopandas Dataframe

texts:
I found this good post on how to get basemap tiles from open sources like OSM, Google, Bing, ESRI using python in geographic lat/long. I'm now trying to add the tiles I successfully got into a matplotlib plot in combination with geopandas dataframe data - which is Projected coordinates.

Getting Geopandas project areas:

aoiD = "o2_boundary.bdy_proj_projectarea" # dataset to query
aoiF = "name"                             # Field to use in query
aoiT = "My Project Site"                  # Value of field I want for query
con = psycopg2.connect(database=databaseCh, user=userP, password=passP,host=hostP)
sql = "SELECT {},geom FROM {}".format(aoiF,aoiD)
gdf = geopandas.GeoDataFrame.from_postgis(sql, con)
minx, miny, maxx, maxy = gdf.geometry.total_bounds

Skipping the part where I convert minx, maxx, miny, maxy to geographic coordinates....but I essentially had to pick the centre of the bounding box but creating an average X and Y in geographic coordinates.

Getting Tiles

a, bbox = getImageCluster(yG, xG, 0.02,  0.05, 14,r"https://mt1.google.com/vt/lyrs=s&x={1}&y={2}&z={0}")
ax.imshow(np.asarray(a),zorder=1,extent=[minx, maxx, miny, maxy])

If I set the extent in imshow to the bounding box of the tiles data, then it does line up, with a bit of error mind. However, I don't want the tile extent - I need to force it to the project area extent. The image should show what happens when I set the extent to the project area min,max X/Y. The imagery is way off.

I couldn't figure out how to change the getImageCluster function to accept projected coordinates. I also can't find any documentation on what the 0.02 and 0.05 figures are for - apart from in the original post "delta_lat,delta_long". Do I need to calculate this delta, the difference of the minx,maxx? I have no idea. The function needs one X, one Y and one Z/Zoom. I would like to provide a bounding box instead. Or request a bigger area, then clip. But I don't know how to clip it and marry it up with the rest of my plot.

Since I couldn't change the request to a bounding box I got different extents of the 2 source datasets :

Geopandas layer /Provided Project Area Box:

X_min: 305293.73475847783. X_max: 316160.562743328. Y_min: 238490.12341525586. Y_max: 245986.05093984905

Bounding Box of Tiled Basemap:

X_min: 309767.457047739. X_max: 314045.453068257. Y_min: 241176.638187177. Y_max: 245650.729011243

And the result is pretty bad, imagery that doesn't line up to the vector data. Any thoughts or different approaches - can the tiled array be filtered with coordinates?


texts: mosaicking and projecting modis land cover data MCD12Q1

texts: I have modis land cover MCD12Q1 data .Its in the sinusoidal projection (.hdf native format ). I have to study the land cover characters from the data for which i have to convert them into .tif format , so that i can use the data in ENVI and arcgis . I used the MRT tool to mosaick and reprojected the tiles into WGS84 Geographic coordinate system converting them into .tif format. When i opened the .tif data , there were just distortion and just black and white strips .I checked the projection type of raw files in arc gis which showed "unkown datum based upon the custom spheroid sinusoidal " . how do i transform the .hdf files to .tif files of valid projection without distortion ?


texts: My QGIS plugin shows GUI but not Works. How to get Error log?

texts:
I am new to PyQgis. I made a simple calculator with PyQt4, now I am trying to run it under QGIS. I connected calculator as plugin, and tested it, its GUI works without error but it does not calculates anything.

How can I get error log?

calculator.py:

# -*- coding: utf-8 -*-

from PyQt4 import QtCore, QtGui
from calculate_ui import Ui_Dialog
from math import *
from qgis.core import *

def derece(degrees):
    return degrees*((2*pi)/400)

class Dialog(QtGui.QDialog, Ui_Dialog):
    def __init__(self, iface):
      QtGui.QDialog.__init__(self)
      self.iface = iface
      self.setupUi(self)
      self.ui = Ui_Dialog()
      self.ui.setupUi(self)
      # Validations
      tempValidator = QtGui.QDoubleValidator()
      tempValidator.setNotation(QtGui.QDoubleValidator.StandardNotation)
      self.ui.at.setValidator(tempValidator)
      self.ui.Rt.setValidator(tempValidator)
      # Signal/slot connections.
      self.setupConnections()

    def calc(self):
      a = float(self.ui.at.text())
      R = float(self.ui.Rt.text())

      T=R*tan(derece(a/2))
      L=(2*pi*R*a)/400
      BS=(R/cos(derece(a/2)))-R

      self.ui.Tt.setText(str(T))
      self.ui.Lt.setText(str(L))
      self.ui.BSt.setText(str(BS))
    def silk(self):
      self.ui.at.clear()
      self.ui.Rt.clear()
      self.ui.Tt.clear()
      self.ui.Lt.clear()
      self.ui.BSt.clear()

    def setupConnections(self):

      self.connect(self.ui.hesaplak, QtCore.SIGNAL('clicked()'),
          self.calc)
      self.connect(self.ui.temizlek, QtCore.SIGNAL('clicked()'),
          self.silk)

You can take a look at the plugin from here.


texts: Converting dxf to shp and keep elevation information QGIS

texts:
I want to create a QGIS plugin that convert dxf files to shp. I didn't use DXF2SHP plugin in QGIS because I have some specific needs. I'm using gdal_ogr library.

dxfdat="*.dxf"
shpLinestring="*.shp"

opt1="-skipfailure  --config DXF_MERGE_BLOCK_GEOMETRIES FALSE --config DXF_INLINE_BLOCKS TRUE -gcp  z  -nlt LINESTRING "
sql1= "-sql \"SELECT *,ogr_style from entities where OGR_GEOMETRY LIKE '%LINESTRING%'\" "
processing.runalg('gdalogr:convertformat',dxfdat , 0,opt1 + sql1, shpLinestring)

I can make the convertion but is it possible to keep the elevation attribute from dxf file ?


texts:
This expression works perfectly when integer fields are 0 instead of NULL but now I am working with a GDB that has many NULL values. I'm not sure how to force CInt to only run if field is not null.

Function FindLabel ( [SQFT_CONC] , [DEPTH_CONC] , [SQFT_ASPHT] , [DEPTH_ASHT]  )
  Dim strInput
  Dim dc
  Dim da
  dc = CInt([DEPTH_CONC])
  da = CInt([DEPTH_ASHT])
  val = dc + da
  if (val = 10) then
     strInput = "IT WORKS"
  end if

  FindLabel = strInput
End Function


texts: What are the main fields of applications of 3D GIS?

texts: I'm doing a scholary research on 3D GIS but i have some problems in finding where they are used. Does anyone know of any real life examples where 3D GIS has been adopted?


texts: WFS Query filter grouping "and" and "or"

texts:
I have a WFS query filtering on a property and a spatial element. Exmaple:

<wfs:Query typeNames="bag:verblijfsobject" srsName="EPSG:28992">
    <fes:Filter>
        <fes:PropertyIsEqualTo matchCase="false">
            <fes:ValueReference>pandidentificatie</fes:ValueReference>
            <fes:Literal>0310100000665851</fes:Literal>
        </fes:PropertyIsEqualTo>
        <fes:And>
            <fes:Intersects>
                <fes:ValueReference>bag:geometrie</fes:ValueReference>
                <gml:Polygon gml:id="geo.1" srsName="EPSG:28992">
                    ...
                </gml:Polygon>
            </fes:Intersects>
        </fes:And>
    </fes:Filter>
</wfs:Query>

Basically this has the structure:

pandidentificatie=0310100000665851
and
bag:geometrie intersects with <gml>

How do I expand with multiple geometries so that I get the structure

pandidentificatie=0310100000665851
and  (
  bag:geometrie intersects with <gml1>
  or
  bag:geometrie intersects with <gml2>
  or
  bag:geometrie intersects with <gml3>
)


texts: Calculate the length of shared boundaries between multiple polygons

texts:
I am trying to calculate the length of shared boundaries between multiple polygons in R. Specifically, I would like to calculate the shared boundary length of all the polygons in Census Place GIS data downloaded from NHGIS from the year 2000. (If you want to obtain the exact data, go to NHGIS.org; Select Data; Under Geographic Levels, select Place; Under years, select 2000; In the table below, select the "GIS Files" tab and download the first file.)

I've looked at this thread, but unfortunately this code throws multiple errors:
https://stackoverflow.com/questions/45338384/calculate-the-length-of-shared-boundaries-between-multiple-polygons. I've tried working with the code to modify it but I haven't been able to get it to work. I believe the problem has to do with the complexity of the intersections. One gIntersection often has multiple classes (Points, Lines, Polygons) with different lengths, as seen in the previous link, but the proposed solution does not work.

My question is identical to this unanswered question:
Length of the boundary intersection between two polygons in R

    require(rgdal)
    require(rgeos)
    library(sp)
    library(raster)
    library(spatialEco)
    library(spdep)

    Shapefile <- readOGR(".","US_place_2000")
    Shapefile <- sp.na.omit(Shapefile, margin=1) # NAs are present for some reason?
    Touching_List <- gTouches(Shapefile, byid = TRUE, returnDense=FALSE)
    # also tried Touching_List <- poly2nb(Shapefile, queen = FALSE)

perimeters <- sp::SpatialLinesLengths(as(Shapefile, "SpatialLines"))

    all.length.list <- lapply(1:length(Touching_List), function(from) {
if(is.null(Touching_List[[from]]==FALSE){
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE)
      if (class(lines) == "SpatialCollections") {
        list.Lines <- lapply(1:length(Touching_List[[from]]), function(to) {
          line.single <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]][to],])
          if (class(line.single) == "SpatialPoints") {
            # Double the point to create a line
            L1 <- rbind(line.single@coords, line.single@coords)
            rownames(L1) <- letters[1:2]
            Sl1 <- Line(L1)
            Lines.single <- Lines(list(Sl1), ID = as.character(to))
          } else if (class(line.single) == "SpatialLines") {
            Lines.single <- line.single@lines[[1]]
            Lines.single@ID <- as.character(to)
          }
          Lines.single
        })
        lines <- SpatialLines(list.Lines)
      }
      l_lines <- sp::SpatialLinesLengths(lines)
      res <- data.frame(origin = from,
                        perimeter = perimeters[from],
                        touching = Touching_List[[from]],
                        t.length = l_lines,
                        t.pc = 100*l_lines/perimeters[from])
      res
    }
})

    all.length.df <- do.call("rbind", all.length.list)


texts: pg_dump Database that was populated with osm2pgsql

texts:
Is it possible to dump a postgres database that was created with osm2pgsql and a .osm.pbf file with pg_dump and import this into a remote postgres database?

After running osm2pgsql the database is 135 Gb.  When I create the dump from postgres, I have an 80Gb file.

When I bzip the dump file I end up with a 17Gb compressed file.

I'm hoping this would be less resource intensive than running osm2pgsql on the remote server.  The remote server has a max of 4 gigs of RAM, and resource limits are easily reached.

The import process took close to 7 hours on a 2.5 ghz MacBook Pro with 16 gigs of ram.  I've been running osm2pgsql for the last two days using --slim and -C 3500.

This is the second time I'm attempting to import on the remote server with osm2pgsql.  I'm thinking that importing a dump of the working database may be more efficient.


texts: Defining custom icon for points in L.geoJSON within fetch call

texts:
I am trying to define a new icon for my points in Leaflet but am drawing a blank. I am using Fetch (see code) to access the data and have it displayed with the default icon but can't work out how to change it to another within the fetch. Thanks.

    let puburl = "http:...
    let pubIcon = new L.icon({
      iconUrl: 'markers/beer.png'
    });

    fetch(puburl)
        .then(function(response) {
            return response.json();
        })
        .then(function(data) {
            L.geoJSON(data, {
                onEachFeature: function(feature, layer) {
                    layer.bindPopup("" + feature.properties.PubName + "")
                }
            }).addTo(map);
        });


texts: ArcGIS OG Cost Matrix not iterating over all origins

texts:
I am working with a large data set of property information (500K+ observations), trying to create a cost matrix to find the drive time of all points to downtown Houston. I have set my origins as the parcel point and have set my destination as a singular point in downtown Houston. I did not change any default settings when I was loading my destinations. When I solve the matrix only 125 observations are returned.

Previous posts about this subject have suggested it's a network connectivity issue. I have tried integrating my network, however that didn't change anything. I have also thought about changing the search tolerance for my origin points, but I don't think that would be a solution.

Large scale example:

Small scale example:

Purple = parcel,
Green point = junction

Looking at the above maps, the analysis is taking place within a small area and avoiding observations which are in the same same neighborhood. Because of these two images I've determined that previous answers to this question may not apply.


texts: How can I avoid incomplete WMS layers when exporting large maps in QGIS 2.6?

texts:
I have a QGIS 2.6 project with around 20 vector layers and an WMS layer (working on Ubuntu 14.04). WMS rendering works OK, although a little bit slow, and exporting A4-size maps also works properly (both to jpg and pdf), but the thing gets tricky when trying to export my final-size map, almost an A0.

Most of the time, some of the tiles of the WMS don't appear in the resulting pdf or jpg (see image below). However, I did export a large a week ago. In fact, the map I cannot export now is exactly the same I printed before, after querying one of the layers.

I have tried several "work-arounds":
 * Print to a pdf printer, instead of exporting it. Same result.
 * Erase the WMS layer and define it again, from scratch.
 * Have deleted the map from the composition and defined it again. No change.
 * Have finally tried with a new project and defining a similar size map...and the problem is still there.
 * I have tried to uninstall and reinstall. Everything seems the same (even the plugins are still there).
 * I have tried to export it in a Windows QGIS (both in a Windows partition and in a Windows virtual machine). The first try with Windows worked OK, but now I face the
   same problem as in Ubuntu.

A couple of ideas:
* could it be a problem of the connection with the WMS server? However, WMS rendering works properly.
* could it be QGIS runs our of memory? If I leave the system monitor open I realize that QGIS would only use 2 of my 8 cores, and not simultaneously, and the core that is in use arrives at 100% of usage. Similarly, used memory never goes beyond 2Gb, but I have 16Gb (and 3Gb in the graphic card).

Hope someone can give me a hint... as I have all my project in QGIS and now cannot get my final result.

Edit
I also tried with limiting tile size, but the problem persisted (with smaller blank tiles, though).
Then I tried doubling the timeout and max-retries parameters, as suggested by . It did work perfectly... at least at that time... One month later I faced the same problem again. This time with a small map (DinA4 size, the previous one was an A0). I have set timeout to 50000 ms and max retries to 20!! No changes.
Two additional details might help figure out where the problem comes from:
* The blank tiles are always the same. Not from a month ago to now, but everytime I face the issue the blank tiles would remain fix during all the trials I do to overcome the problem. That doesn't point to a problem with the internet conection, does it?
* The export image is created in a few seconds, not more than 30, which also suggests timeout parameter is not involved (or not working).

Edit
I found that when this problem arises it stays there no matter what I do, with always the same tiles missing. However, if I let it for a couple of days, it disappears. I mean, once it happens it will happen again in all trials during a couple of days; but if I forgot about it and try two days later it works perfectly again. Seems as if it has some kind or memory... Once the communication has failed for whatever reason it will remember it for a couple of days...


texts: Interchange y, x to x, y with geopandas (python) or QGIS

texts:
I realized that the order of my geopandas dataframe is wrong, it has y, x instead of x, y.

This makes it difficult to use with other modules that needs exactly x, y type of coordinates. Is there any function in geopandas to interchange the coordinates without the need to enter in each geometry element, geometry type and change coordinates one by one? If not is there a shapely function to do it similarly instead of collecting all points one by one and change them?

More info: I have a geodataframe composed of a lineStrings (shapely object) in the geometry field for each row. If it's easier to do it in QGIS, I could also do it there.


texts: Using compiled python file for QGIS script PyQGIS

texts:
Is possible to use compiled python file (*.pyc) for "Open Existing Script ..." in the Scripts Processing module?
File dialog does not enable selecting *.pyc files.
I want to use QGIS scripts as PYC files.


texts:
I'm fairly new to working with versioned geodatabases and have found that one of my versions has over 600 conflicts which require review and replacing.

I suspect that all objects will need to be replaced with the Conflict Version, but I would like to confirm my suspicions.

Is there a quick / automated way to review all of them without going through one by one?

I'm relatively new to learning Python too.


texts: How to show what enable_outdb_rasters is set to in Postgis

texts:
I know that you can enable out of database raster support via the query:

SET postgis.enable_outdb_rasters TO True;

Before running this query I'd like to verify that postgis.enable_outdb_rasters is set to False. Is there a query to show or select what the current postgis.enable_outdb_rasters variable is set to?


texts:
I'm trying to do some fire spread simulations in QGIS with GRASS fire spreading simulation module. QGIS ('Processing Toolbox' > 'Geoalgorithms' > 'Misc' > 'Fire spreading simulation') GUI.
Already searched in help files (insufficient information), asked, googled...

Is there anyone that knows how to work with this module?
Can you explain (or point to a tutorial on) how to run the module and what are the input formats?
(ASCII, GeoTiFF, can it be .shp for ignition?, fuel model defs: is it FARSITE/FlamMap fuel models or personal fuel models?)


texts: Can you edit a QGIS layer style file (.qml) to create random colors?

texts: I have a .qml file that I am using on a project importing photos using the 'import photos' plugin.  Every import batch needs to have the color manually changed after importing in order to distinguish it from the other imports.  This is something I'll be doing hundreds of times.  Is there a way to set the qml file to randomize the color when importing the data?


texts: Reading zipped gpkg files (i.e. `*.gpkg.zip`) with Python

texts:
I have many large gpkg.zip files that I would like to read directly into Python (to then upload to BigQuery). I can unzip each, read with geopandas, then upload -- however, since unzipping triples the space the file takes up, it starts to really clog up my computer's relatively small hard drive.
Is there a way to read a compressed file in directly?


texts: Adding a tab to the QGIS project properties

texts:
I would like to know, if there is the possibility to add new tabs to the QGIS project properties? (for example, under the "QGIS Server")

The goal is to have the advanced options of a QGIS plugin directly in the project properties.


texts:
I have several raster files that are basically different indices produced with the same DEM. I want to visualize the raster layers in a way that the audience can see that they are stacked.

I am wondering is there a software that produces this effect?.

An example is given below (From Neon).


texts: Represent joined data in Info window (ArcGIS JSAPI 3.1)

texts:
Currently I am facing following problem

I have joined  the SQL Server data base to my shape files (In ArcGIS Desktop using “join” option)and published on ArcGIS Server and working on REST service in JavaScript API.
After joining the tables to the map… all columns are changed. Let us say info.id, attr_name.ratio etc.

Therefore, when I am using info window option I cannot use “info.id” it is considering only “info” it gives an error if I used “info.id”.

Any help will be great...

Software’s: ArcGIS DT 10.1 and Server 10.1

Code Sample:

'var featureLayer = new esri.layers.FeatureLayer("http://<myserver>:6080/arcgis/rest/services/type/MapServer/1",

{
     // mode: esri.layers.FeatureLayer.MODE_SNAPSHOT,
      mode: esri.layers.FeatureLayer.MODE_ONDEMAND,
      infoTemplate:popupTemplate,
      outFields: ['info','attr_name']

}); '


texts:
I would like to add halos on my labels on an Openlayers map.

I have tried adding the following to my code but no halos display.

labelHaloColor: "#FFFFFF",
labelHaloWidth: "2",

I got that from here:

http://trac.osgeo.org/openlayers/ticket/2965

Do I need to use a newer version of OpenLayers to get this working? I am currently using the dev version.


texts: ogr2ogr specify lat/long order on conversion to GeoJSON

texts:
When using ogr2ogr to convert a polyon from a .shp file to GeoJSON, I noticed that it writes the coordinates in order lat, long. I noticed this because whenever I tried to view the GeoJSON on a map, it would place the polygon in the middle of nowhere; I then read a bit more about GeoJSON and realized that the format requires the (x,y,z) order (long,lat,altitude) [source]. Is there a way for me to specify that ogr2ogr should write the GeoJSON in that order? Or is there any other way for me to get around this problem?

I'm 100% open to the possibility that I'm misunderstanding GeoJSON and that it doesn't actually have to be in long,lat,altitude order.


texts: I have a polyline layer has a unique text field that identifies a path.   I have made a union with a separate 1km lattice grid.  I am needing a way to get a "count" of the number of uniquely identified paths that fall within each 1km grid spacing.  Is there a tool in ArcGIS (of note I have an ArcInfo license) to do this?


texts: How to show raster's cell value over its map in levelplot?

texts:
Hy, Is there any way to plot each cell of raster with its values in level plot ?
Raster is attched here(
)
I have tried these https://stackoverflow.com/questions/6569666/showing-data-values-on-levelplot-in-r, but it is only ploting the matrix values.
If anybody can help in this regard.


texts: How can I eliminate leftover border/boundary in QGIS?

texts:
I am trying to generate a shapefile which contains the unincorporated areas of St. Louis County, Missouri.

I'm using two shapefiles: St. Louis County from the 2012 census, followed by St. Louis County places from the 2012 census. I have removed all CDP from the places file, so that it consists only of incorporated municipalities.

Here is St. Louis County from the US Counties shapefile:

Here are the incorporated municipalities from the Places shapefile:

In QGIS, I dissolve the municipalities together to make one "incorporated" polygon:

Finally, in QGIS I take the difference of incorporated polygon from the St. Louis County polygon. In theory, this should yield my unincorporated St. Louis County polygon.

As you can see, it's almost there.

Apparently the outer borders of the Places shapefile and the St. Louis County shapefile were slightly off, so that when I took the difference, some tiny slivers of the outer edge remains.

And so my unincorprated shapefile ends up with a sort of leftover outline around empty space, where there should be no border at all.

How can I eliminate this extra border/boundary/leftover thing in QGIS? Can anyone help?


texts: QGIS - Adding Raster Data to a Scanned Map

texts:
I'm a bit of a novice to QGIS, so apologies if this is a silly question:

I've scanned a geological map into QGIS, which shows different types of rock in different colours. It's a raster layer, and I'd like to have it work like the elevation layer I've imported from the USGS - that is, for each pixel to have a value (for example, 'chalk') so that I can then tabulate the type of rock on which various features are situated. How do I go about this?


texts: Very large VRT file generated by gdal after cutline

texts:
In Python, I generate a virtual raster using gdal warp with a cutline as per below

gdal.Warp(warpedFilePath, ds, format = 'VRT', cutlineDSName = csvFilePath, srcNodata = 0, dstAlpha = True, cropToCutline = True, dstSRS = 'EPSG:4326')

If the cutline is within the image, the vrt file is of "normal" size (a few KB) and the processing afterwards (when turning it into a GTIFF) is fast.
However, if the cutline has matching edges with the image then I get a huge VRT file, eg several MB and the processing is over an hour. Looking at the Cutline value, I can see that every border pixel is being added! when the translate is done to have a tiff file, the process takes over an hour...

Is this a "normal" behaviour? is there an option to reduce the number of points in the cutline?

My gdal version is 2.2.4
The input raster is in UTM and the cutline is a csv file in 4326.

Edit:
CSV file for the cutline

    ID, WKT
1, "POLYGON ((-4.68694239006474 54.9863989533752,-4.68563773368381 54.9553739114236,-4.92089428671608 54.9518432821996,-4.948233 55.021719,-4.711947 55.051303,-4.68694239006474 54.9863989533752))"

Unfortunately I cannot provide the image, for one it is too large, around 1GB, and it is bound to its commercial terms. That said I am quite convinced the same would happen with a free Landsat image for instance.


texts: Sort and subset feature classes using Python script

texts:
Here's my situation: I've got multiple utility feature classes, and each fc is subtyped according to one the primary attributes, which is a String. Each fc has at least 12 attribute fields, and each attribute is formatted as "String". My goal is to write a Python script that will sort, separate and export out feature classes for each subtype of the data. Does anyone have a good place to start or an idea of how to set-up the skeleton of the Py to get my script going?
As an example (or reference) I'd like to pull all,  "Watermain, 8in, PVC, Owner" out of a jumble of utility data.


texts: Correctly Setting Spatial Reference for Geometry Objects when Geoprocessing

texts:
I'm attempting to find the nearest distance between a set of points, and a set of polygons. I don't want to use near tables because I have to do this operation for a large number of different selections on my polygon file, and I don't want to generate that many tables - all I need is a number. I'm trying to do this using geometry objects and the distanceTo() method, and storing a variable that keeps the smallest value for distanceTo(). In the end, it should just give me a single number: the shortest distance.

My problem so far has to do with spatial references I think. My point file is in UTM, and my polygon file is in state plane. Originally, the numbers I was getting were too large - millions of feet apart. So, I grabbed the spatial reference from my point file, and inserted it into both of my geometry objects when using the searchCursor. But now the numbers being returned are all 0 - and I know that my points don't all intersect.

I also know that all my selections on the polygon file are working properly - I've checked them ad nauseum - so I really don't think I'm getting erroneous numbers because of that.

So, my question is, have I entered the spatial references correctly? And if so, what else could be causing this issue? Is there something (there are probably many things...) that I don't understand about the distanceTo() tool?

shortestDistance = 10000000

pointFile = r"c:\data\points"
desc = arcpy.Describe(pointFile)
sr = desc.spatialReference

polyFile = r"c:\data\polys"
#these files get turned into feature layers
#in order to perform spatial selections
#and attribute queries

for row in arcpy.da.SearchCursor(polyFile, ["SHAPE@"], spatial_reference=sr):
    for part in row[0]:
        for pnt in part:
            if pnt:
                x = pnt.X
                y = pnt.Y
                point = arcpy.Point(x, y)
                array.add(point)
            else:
                arcpy.AddMessage("Polygon with Interior Ring")

    vegPoly = arcpy.Polygon(array)

for row in arcpy.da.SearchCursor(pointFile, ["SHAPE@XY"], spatial_reference=sr):
    #make it a point geom
    x, y = row[0]
    point = arcpy.Point(x, y)
    centroid = arcpy.PointGeometry(point)

    distance = centroid.distanceTo(vegPoly)

    if distance < shortestDistance:
        shortestDistance = distance


texts: How to enable the snapping options for vector layer using the python programming?

texts: I  want to enable the snapping options for the the existing layers using python programming. How can I do this?


texts:
I found this code from Google Earth Engine.
/*
Author: Sofia Ermida (sofia.ermida@ipma.pt; @ermida_sofia)

this function mask clouds and cloud shadow using the Quality band

to call this function use:

var cloudmask = require('users/sofiaermida/landsat_smw_lst:modules/cloudmask.js')
var TOAImageMasked = cloudmask.toa(image)
var SRImageMasked = cloudmask.sr(image)
or
var TOAcollectionMasked = ImageCollection.map(cloudmask.toa)
var SRcollectionMasked = ImageCollection.map(cloudmask.sr)

INPUTS:
        - image: <ee.Image>
                image for which clouds are masked
OUTPUTS:
        - <ee.Image>
          the input image with updated mask
*/

// cloudmask for TOA data
exports.toa = function(image) {
  var qa = image.select('BQA');
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image.updateMask(mask);
};

// cloudmask for SR data
exports.sr = function(image) {
  var qa = image.select('pixel_qa');
  var mask = qa.bitwiseAnd(1 << 3)
    .or(qa.bitwiseAnd(1 << 5))
  return image.updateMask(mask.not());
};

I am interested in using it through the Python API, but I have a problem in the understanding of the syntaxis. I am totally novice working on JavaScript - Node JS and I have issues on how to implement this kind of statement: exports.sr = function(image) in Python. What would be the equivalence in the Python API?


texts: ArcGIS API for Silverlight : query that displays layer attributes in a dataGrid

texts:
I'm working with ArcGIS API for Silverlight to build a webmapping application and I'm trying to display my attributes in a datagrid. I made a try by using this code :

void QueryTask_ExecuteCompleted(object sender,
         ESRI.ArcGIS.Client.Tasks.QueryEventArgs args)
  {
   FeatureSet featureSet = args.FeatureSet;

  if (featureSet != null && featureSet.Features.Count > 0)
        QueryDetailsDataGrid.ItemsSource = featureSet.Features;
  else
        MessageBox.Show("No features returned from query");
  }

  private void QueryTask_Failed(object sender, TaskFailedEventArgs args)
  {
        MessageBox.Show("Query execute error: " + args.Error);
  }

But I get the result No features returned from query. I'm not sure, but it seems like the condition featureSet != null && featureSet.Features.Count > 0  isn't valid.


texts:
I have one shapefile which is divided into multiple villages (polygons). I can see the name of the villages from attribute table of the shapefile. For each village, I want to extract elevation from ASTER DEM image. In other words, I want to extract based on the attribute (ObjectID/Name) of one shapefile.

Is there any efficient way (i.e. coding, model building) to do it? [I am using ArcGIS 10.5.1]


texts: SAGA raster calculator, ifelse command syntax problem


texts: Do postgres 9.1, postgis 2.0 from stackbuilder and pgrouting 1.0.3. work together?

texts:
I am new to pgrouting and I would like to check if the version I have installed works fine with postgres and postgis.

I have installed postgres 9.1, postgis 2.0 from stackbuilder and pgrouting 1.03.

Any suggestion how I can check if pgrouting works fine?


texts: Converting from custom projection northing/easting coords to lat/long in decimal degrees

texts:
I've got a custom projection file that's applied to all data (per request of a client), and another team is requesting ids/coordinate pairs for assets. I tried ripping out X/Y off the point data, but I noticed the coordinates make no sense - they are not standard UTM northings/eastings.

For instance:

A point with real-world coords of 43.814794N, -76.02576907720452 has these coordinates in the weird custom projection: X: -534643.85116969, Y: 303444.080342568

The GCS is WGS84, but none of the online converters or formulas I can find work, and nobody has any clear instructions on how to create a formula. What I need is to write a Python function that will convert the eastings/northings to decimal degrees. Anyone have knowledge of how to develop such a formula?

Here's the info on the projection:

  False_Easting       |  0.000000000000000000

  False_Northing      |  0.000000000000000000

  Central_Meridian    |  -74.000000000000000000

  Standard_Parallel_1 | 42.000000000000000000

  Standard_Parallel_2 | 44.000000000000000000

  Latitude_Of_Origin  | 43.000000000000000000

And here's the info on the GCS (WGS84):

Name: GCS_WGS_1984

Angular Unit: Degree (0.0174532925199433)

Prime Meridian: Greenwich (0.0)

Datum: D_WGS_1984

  Spheroid: WGS_1984

    Semimajor Axis: 6378137.0

    Semiminor Axis: 6356752.314245179

    Inverse Flattening: 298.257223563


texts: Adding picture to ArcMap PageLayout using ArcObjects in .NET?

texts:
I'm using ArcGIS Desktop 10.6 and I'm creating an Add-in to simultaneously add text (textelement populated from user input) and an image (company logo).  The issue i'm having is with the image.  I know there is a pictureelement but I have no idea how to work with a stdole.IPictureDisperse so I'm trying to see if I can use a PictureMarkerSymbol but no luck.

System.Uri myPictureUri = new System.Uri(picUrl);
PictureMarkerSymbol picm = myPictureUri as PictureMarkerSymbol;
graphContainer.AddElement(picm as IElement, 0);


texts:
I use the following code to calculate the areas of a multi-polygon shapefile. However, to validate these results, I calculate the polygon’ areas directly in table attribute by creating a new field and applying the formula: $area. I found discrepancies in calculated areas by both methods.

Does the QgsDistanceArea class consider the layer CRS (in my case I use EPSG:32632) and $area does not, that means the geographic WSG84 EPSG:4326 is considered by default?

def area_poly(layer):
  L=[]
  features=layer.getFeatures()
  d = QgsDistanceArea()
  d.setSourceCrs(layer.crs(), QgsProject.instance().transformContext())
  d.setEllipsoid(QgsProject.instance().ellipsoid())
  for i, element in enumerate(features):
    for poly in element.geometry().asMultiPolygon():
        m = d.measurePolygon(poly[0])
        L.append((i, round(m,1))
  return L


texts: Shapes don't match google maps QGIS 2.18.7

texts:
I have some shapes (in EPSG:4326 wgs 84) that match with Openlayers ->google street maps only when I do enough zoom-in.

But When I zoom-out seems that my shapes move to the east.

What could be happening?

I think it might be something related to projection. But since I am newbie in GIS I dont know deeply qgis tools and configurations.

In picture below Black line should be matching with red line, than happens only if I zoom-in.

By the way, the more zoom-out, the more shapes move to east.


texts: Calculating flow distance of pixel from outlet using R?

texts:
I have created these raster files for a river basin:

DEM
Flow Direction
Flow Accumulation

Also, I have another raster with values of rainfall.

All the 4 rasters have same spatial extent, resolution and dimensions.

How can I calculate the flow distance of each pixel of rainfall from the outlet of the basin?

Theoretically, I have to follow the Flow Accumulation values, but I'm not sure how to implement it and whether there are existing implementations.

I am currently working with raster package in R, but I am open to other platforms.


texts: Help with Point on Polygon overlay in ArcPy

texts:
I am trying to do this simple task but not sure how to do this in Python. I have a point shapefile and I want to overlay with polygon layer to get all the polygons for each point. After this would like to loop through all the rectrods to get the attributes. Appreciate any feedback.

wellRows = arcpy.UpdateCursor("H:\points.shp")

for wellRow in wellRows:
    wellPoint = wellRow.SHAPE #Not sure how to have shape as input...
    landRows = arcpy.SearchCursor("H:\Poly.shp", "", "", "","")

    for landRow in landRows:
        print landRow.NAME


texts: Can anyone figure out this projection? Eckert III looks close but is not quite right.


texts: Check whether a point is inside a polygon or not in MySQL

texts:
I am trying to get the points those are within the shapefile. I have tried the Spatial functions and it's showing an error:

3033 - Binary geometry function st_contains given two geometries of different srids: 4326 and 0, which should have been identical.

My query is
SELECT
lightning.lightning_id, lightning.latitude, lightning.longitude,
lightning.flash_type,
DATE_FORMAT(lightning.lightning_time, '%D %b %Y %T %p') as lightning_time
FROM lightning, districts_s
WHERE
ST_CONTAINS(districts_s.SHAPE, Point(lightning.longitude, lightning.latitude))
AND lightning.height > 50

To rectify the error I got a solution from here and this link but no result, as this also showing an error:

ST_Transform does not exist

update 1
ERROR ON ST_SetSRID


texts: Arcobjects Multithreading in C++ - 'Aborted' when modifying Geometry

texts:
My project is required to be cross-platform, so I am using the C++ API for ArcObjects. My application creates a new thread which listens to a network socket broadcasting location information. I am attempting to read in the information and add points to my DynamicLayer.
My problem comes when I try to use a new IPoint in the second thread. I can make the call:

IPointPtr newpoint(CLSID_Point);

But on the next line when I execute:

newpoint->PutCoords(x,y);

The application fails with a nonspecific error 'aborted'.

Relavant information: I am creating a new thread by calling pthread_create on a static method.

I have read the ESRI document on multithreading (which like the rest of the documentation is woefully absent of c++ specific information), but didn't see anything I thought could be the problem. Any thoughts?


texts: Why might my MapServer WMS layer be invisible in QGIS if they work fine in a browser?

texts:
I've managed to create a WMS feed from Mapserver showing a tile index to some rasters and a shapefile and it appears in my browser OK. If I go to QGIS (1.8.0 Windows StandAlone) and I connect to my WMS feed I can see the layer names appear in the dialog box. However when I add the layers to the map I don't see any layers. I tried zooming to them etc but nothing is visible on the screen.

Does anyone have any ideas why this might be?


texts: Processing raster present as MS Access database

texts:
I wanted to know if anyone has experience working with WISE soil database, I am trying to process world soil data which is downloaded from here

https://data.isric.org/geonetwork/srv/eng/catalog.search#/metadata/dc7b283a-8f19-45e1-aaed-e9bd515119bc

After I downloaded the folder, there are 5 files present:

hwsd.bil
hwsd.blw
hwsd.hdr
HWSD (Microsoft Access database)
HWSD_META (Microsoft Access database)

I have not really worked with raster data present an an MS Access database. What I need to do is for a given lat long, extract all the data from the above soil database. How do I go about doing this in R? Is there an online tutorial?


texts: gdal_sieve in QGIS Python script not running?

texts:
I have a raster that I want to filter using gdal_sieve.py. However I'm having some problems with my syntax on the last line.

from qgis.core import * #importing the QGIS processing core
from qgis.utils import iface #importing the iface
import os, subprocess

subprocess.call(["gdalsieve.py", "-st", "300", "8", \
r"/data/final_classification.tif", \
"-of", "GTiff", r"/data/brent/LUCID/sieved.tif")]

  File "input", line 3
      "-of", "GTiff" r"/data/brent/LUCID/buk.tif")] ^ SyntaxError: invalid syntax


texts:
I need to aggregate disjoint polygons within a certain distance in ArcGIS but it seems that the perfect command to do it, Aggregate Polygons, is only available with an ArcInfo that I don't have.

Would you know whether there's another way to do this with another command?

Below is a sample of two polygons to be merged.

NOTES:

the usual method using Editor will not work because it's many polygons
Resampling using Majority will not only do either because it increases the area of the polygon too much.


texts: How to configure osm2po to use the specified maxCost in TspDefaultMatrix

texts:
The problem

I am using osm2po's TspDefaultMatrix to create a cost (time) matrix:

new TspDefaultMatrix(graph, vertexIds, maxCost, Log.stdout(Log.LEVEL_LOG), properties)

Where properties contains:

findShortestPath: false
ignoreRestrictions: false
ignoreOneWays: false
heuristicFactor: 1.0

I tried many value for maxCost (0, 1, 2, Float.MAX_VALUE). Except for maxCost = 0, I always get this in the log:

DEBUG  Heuristic MaxCost set to 0.2134047

followed by some:

WARN   Matrix [8][9] unreachable!
WARN   Matrix [8][3] unreachable!

because, indeed, some of the vertexId are more than 0.2134047 hours appart.

The question

I want to know how this value (0.2134047) is calculated and how I can effectively tell TspDefaultMatrix to use my maxCost.

By the ways, I did not found (was not searching a lot) the osm2po source code. Is it available ? I would have been able to answer my question myself.

EDIT:

By the way, I was using version 4.8.8 of osm2po to build the .gph file and at runtime.


texts:
I'm very new to ArcGIS having used Global Mapper in the past to do simple viewshed analyses of factory stacks at a given radius.

Now I need to produce a 2km radius ZTV of two buildings next to each other each of a given height.

I have the footprints in CAD saved out as DXF (though I realise this may be bad practice?!) and the DTM in place (UK Ordnance survey ASC tiles), but I cannot work out how to extrude those buildings in ArcGIS and perform a viewshed analysis on them?


texts:
I have Android/Java code that adds layer to a MapView, calling a RESTful service, like so:

mDynamicServiceLayer = new ArcGISDynamicMapServiceLayer(getResources().getString(R.string.some_url));
mMapView.addLayer(mDynamicServiceLayer);

This works fine. I want to also add a layer that is an mxd file on the device

<string name="mxd_file">file://storage/emulated/0/Download/mymap.mxd</string>

mDynamicServiceLayer = new ArcGISDynamicMapServiceLayer(getResources().getString(R.string.mxd_file));
mMapView.addLayer(mDynamicServiceLayer);

This throws the exception:

java.lang.IllegalStateException: Scheme 'file' not registered.

How do I do this?


texts:
I am trying to rotate a Polygon type(see below) with the center of rotation(not necessary in the center) but somewhere within the boundary. How can I do this?

I am looking at rotate and affine_transform. The polygon vertex point (4 corner) is in meters (dimension) which i will need to transform it to lat lon in EPSG:4326 later stage

I was tempted to use rotate since parameter to use as it is simple ie rotate(geom, degree, centerofrotation) but was told that  it does not work for lat lon coordinates.

To use affine transform, can anyone show me how to calculate for the transform matrix?

for example, the polygon is

POLYGON ((11576467.24265137 150763.0256523852, 11576467.24265137 150805.3789849052, 11576509.04265061 150805.3789849052, 11576509.04265061 150763.0256523852, 11576467.24265137 150763.0256523852))

and center of rotation (11576488.142650992, 150779.6256524612)

Thanks


texts: Display_expression function on QGIS 3.14

texts:
I would like to know what exactly do the function display_expression introduced in QGIS since version 3.14.
The 3.14 Changelog is not clear on it and also reading the GitHub commit didn't help me to solve my doubts about the use of this function.


texts: Re-ordering bar chart from Data Plotly in QGIS 3.14.16

texts:
I have successfully created a bar chart using Data Plotly within QGIS 3.14.16. The x axis consists of 17 place names (NAME) and the Y axis shows the value (VALUE IN £P) of one particular lease at each location, ranging from £0 to £46. This is an historical chart. However I would like the plots to be ordered by value so that the lowest value is on the left and the rest graduated to the highest value on the far right of the plot.
Is this possible? Can any kind person provide me with an expression to achieve this please?


texts: Why QgsMapCanvas only appear after lost focus?

texts:
I created the simple frame .UI file, compiled with pyuic4, and tried the ShapeViewer.py to load the shapefile "world_borders.shp". However, the map doesn't display at the first time. I have to click out the window to make it appear. After everything works!
I am working in Windows 7.
Instalation "osgeo4w-setup-x86".
Qgis 2.2.0-Valmiera

Why does QgsMapCanvas in Qgis 2.2.0-Valmiera only appear after losing focus?

Here is my code:

class ShapeViewer(QMainWindow, Ui_MainWindow):
    def __init__(self):
        QMainWindow.__init__(self)
        self.setupUi(self)
        self.canvas = QgsMapCanvas()
        self.canvas.useImageToRender(False)
        self.layout = QVBoxLayout(self.frame)
        self.layout.addWidget(self.canvas)
        file = QFileDialog.getOpenFileName(self,"Open", ".", "Shp(*.shp)")
        layer = QgsVectorLayer(file, "My Shape", "ogr")
        QgsMapLayerRegistry.instance().addMapLayer(layer);
        self.canvas.setExtent(layer.extent())
        cl = QgsMapCanvasLayer(layer)
        self.canvas.setLayerSet([cl])
        self.show()

def main(argv):
    # create Qt application
    app = QApplication(argv)

    # Initialize qgis libraries
    QgsApplication.setPrefixPath(qgis_prefix, True)
    QgsApplication.initQgis()
    # create main window
    wnd = ShapeViewer()
    wnd.show()
    # run!
    retval = app.exec_()

    # exit
    QgsApplication.exitQgis()
    sys.exit(retval)
if __name__ == "__main__":
    main(sys.argv)


texts: How to convert a bing maps quadkey to mercator coordinates

texts:
How do you convert a bing maps quad key, such as "03200320023" to coordinates in mercator virtual earth (epsg 3857) coordinate system?

I've got an algorithm to convert to coordinates in latitude longitude (e.g. the quad key above comes to "-84.19921875,33.7243385314941,-84.0234375,33.8704147338867"). I could then convert from latitude longitude to mercator coordinate, but that smells far too expensive in processing power.

The mercator virtual earth result I'm looking for given the above quadkey is: "-9373014,4011415,-9353446,4030983"

Also a question on stackoverflow


texts:
I use Orfeo ToolBox(OTB) in conjunction with QGIS 3.4.
So. I use the qgis-otb plugin.

I was not a problem when I linked OTB with my company's PC,
but when I tried the same way on my personal laptop, I got the following error.

  AttributeError: 'NoneType' object has no attribute 'setProgress'

Traceback (most recent call last):

File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\gui\ConfigDialog.py", line 86, in apply
    self.config_widget.accept()
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\gui\ConfigDialog.py", line 325, in accept
    setting.setValue(str(self.items[setting].text()))
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins\processing\core\ProcessingConfig.py", line 293, in setValue
    self.validator(value)
  File "C:\qgis-plugins\qgis-otb-plugin\otb\OtbAlgorithmProvider.py", line 267, in validateAppFolders
    OtbUtils.executeOtb(commands, feedback=None)
  File "C:\qgis-plugins\qgis-otb-plugin\otb\OtbUtils.py", line 76, in executeOtb
    feedback.setProgress(0)
AttributeError: 'NoneType' object has no attribute 'setProgress'

Is there any solution?


texts: What does the "trusted" flag mean for QGIS plugins?

texts:
Can anyone provide some insight into the meaning of the "trusted" flag (true/false) that is listed on the QGIS plugin website when viewed here -

http://plugins.qgis.org/plugins/plugins.xml?qgis=2.14

Our IT department have stumbled across this page rather than the usual
webpage view (http://planet.qgis.org/plugins/) and now they are concerned that those flagged as "Trusted = false" are potentially dodgy.

We found this page that provides some explanation - http://hub.qgis.org/projects/quantum-gis/wiki/PluginRepository

But can someone please clarify if a FALSE trusted rating simply means the plugin needs to be reviewed and approved by moderators before it will be made visible to the public or if it actually means there is a level of uncertainty even once it is approved and appears on the public plugin lists?


texts: Manually classifying graduated data outside of range in ArcGIS for Desktop?

texts:
I have two shapefiles representing grids with values that are symbolized with graduated colors, shown below:

You can see that the two shapefiles encompass similar ranges of data, but assign slightly different colors to ranges in these data.  I would like to symbolize the ranges in each of the two shapefiles using the same colors for the same ranges, enabling easier comparison between the files and allowing a single legend to be used.

However, when I attempt to classify the ranges with a manual interval, ArcGIS forces what should be the uppermost range to also include the minimum value from the dataset.  You can see this in the "-81.64 - 10.00" range of the left data frame, below.  This forces the entire layer to be symbolized with this color.  Essentially, ArcGIS is requiring that the minimum value in the dataset be used as a value in one of the ranges.

Is there a way to circumvent this feature?


texts: Is possible to create a shapefile with a specific encoding using ArcObjects?

texts:
I have the code to create a shapefile programmatically, something like this 28310 - Create a new empty shapefile programmatically

moreover I am able to put the data from a layer in it.

It works perfectly.

Some client hands to me several layers that have data in Russian. When I export from my application to a shapefile I got the same behavior described here esri - Character encoding support in geodatabases and shapefiles - Geographic Information Systems Stack Exchange .

the fields with Cyrillic alphabet appears like ???????? .

I found this article: 21106 - Read and write shapefile and dBASE files encoded in various code pages  but it does not apply to ArcEngine (ArcObjects) .

Is there a way to change the encoding when I'm creating the shapefile using the workspace provided by the shapefileworkspacefactory?


texts: Landsat 8 cloud/cloud shadow/snow masking in Google Earth Engine

texts:
I am a new Google Earth Engine user and I am trying to understand basic procedures/functions. More specifically, I want to remove clouds, cloud shadows and snow pixels form a Landsat 8 collection and then, calculate the NDVI. I have found the following sample code but I can not understand how to remove apart from the clouds, the cloud shadows and the snow. I think the Fmask might be he way to do this but I can not understand this.

Example code from Google Earth Engine API:

// This function masks clouds in Landsat 8 imagery.
var maskClouds = function(image) {
  var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
  return image.updateMask(scored.select(['cloud']).lt(20));
};

// This function masks clouds and adds quality bands to Landsat 8 images.
var addQualityBands = function(image) {
  return maskClouds(image)
    // NDVI
    .addBands(image.normalizedDifference(['B5', 'B4']))
    // time in days
    .addBands(image.metadata('system:time_start'));
};

// Load a 2014 Landsat 8 ImageCollection.
// Map the cloud masking and quality band function over the collection.
var collection = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
  .filterDate('2014-06-01', '2014-12-31')
  .map(addQualityBands);

// Create a cloud-free, most recent value composite.
var recentValueComposite = collection.qualityMosaic('system:time_start');

// Create a greenest pixel composite.
var greenestPixelComposite = collection.qualityMosaic('nd');

// Display the results.
Map.setCenter(-122.374, 37.8239, 12); // San Francisco Bay
var vizParams = {bands: ['B5', 'B4', 'B3'], min: 0, max: 0.4};
Map.addLayer(recentValueComposite, vizParams, 'recent value composite');
Map.addLayer(greenestPixelComposite, vizParams, 'greenest pixel composite');


texts: v.dissolve (and other scripts) not in GRASS path

texts:
I see that I have the v.dissolve module in my grass64/scripts dir, but it is not found in my path from the GRASS command line.  I have the GRASS_ADDON_PATH environment variable set:
GRASS_ADDON_PATH=/usr/local/grass-6.4.1/scripts/

And the module itself is located in:
/usr/local/grass-6.4.1/scripts/v.dissolve

What am I missing?

Thanks,

Roger


texts: ArcGIS JavaScript API: How to add a CSV as a Layer?

texts:
I've created a map on ArcGIS online, adding a layer by uploading a csv file. This works fine, see map at link http://www.arcgis.com/home/webmap/vi...1ba43a0256c06b.

Now I'm trying to simulate this in JavaScript API. First I create a Web API service that returns the same data set I've got in csv file.
 This is published on my local webserer. In my web page, I create a webmap, adding a layer that references this web API:

 function init() {
 var options = {
 basemap: "gray",
 center: [153.115672, -30.296312],
 zoom: 14
 }
 var mapServiceURL = "http://[myserver]/TRWDataManagerAPI/API/OzRegs";

 // Create map and call Init function
 map = new esri.Map("mapDiv", options);

 map.addLayer(new esri.layers.ArcGISDynamicMapServiceLayer(mapServiceURL));
 }
 dojo.ready(init);

The webmap is displayed, but the layer doesn't have any effect. What's missing? Am I at the right track? Your help is greatly appreciated.


texts: Connect pg routing to HTML webpage for routing application

texts:
I am trying to create a web based routing application for my final year project. I am currently trying to connect my web page with pg_routing. I have all my features loaded in GeoServer and PostgreSQL. I can connect my features to my web page using WMS from GeoServer (example below), but pg_routing does not function in the same way.

Is there a way to do it and if so, how?

var Buildings = L.tileLayer.wms('http://localhost:8010/geoserver/wms' , {
    layers: 'IST_Mosaic:Buildings',
    format: 'image/png',
    transparent: true});


texts: I want to cut a polygon layer with many small polygons (green) to the borders of a polygon-layer with only one large polygon (red). I use the standard geoprocessing clip tool in ArcMap 10.4. Why is it not clipping to the borders?


texts: Saving Coordinates and Depth in a File on click

texts: I am working on Hydrographic maps and I have used the interpolation plugin and got the depth at all points in a raster layer. Now I need to save the Lat\Long and depth in a text file. I am aware of the copy_coords plugin but it saves only x and y coordinates. Similar to copy_coords, I need a code/plugin which also copies (or saves) the depth.


texts: Calculating loss area in a specific elevation class in Google Earth Engine

texts:
I would like to calculate the loss area in a certain elevation class in GEE. When I run the code 1 below, it gives the same amount as the total loss area in my study region, code 2. I switched arealoss and class3 in code 1 as well, and didn't work. How can I calculate the loss area in a certain elevation class?
code 1:
var arealoss= loss.multiply(ee.Image.pixelArea()).divide(10000);
var class3 = elevation.eq(3).selfMask();
var stats1 = arealoss.reduceRegion({ reducer:
ee.Reducer.sum(),
geometry: class3.geometry(),
scale: 30,
maxPixels: 1e9,
bestEffort: true });

code 2:
var stats2 = arealoss.reduceRegion({
            reducer: ee.Reducer.sum(),
            geometry: peru.geometry(),
            scale: 30,
            maxPixels: 1e9,
            bestEffort: true });

Besides, I want to repeat this calculation for 7 different elevation classes. is it possible to write a function for this calculation in GEE?


texts: Best open source operating system for open source GIS?

texts: I'm wanting to get into open source GIS, and it looks like PostGIS on Postgres is a good DBMS to go with, and there are plenty of posts on desktop editors. What I haven't come across (and forgive me if I just didn't find it) is info on what open source operating systems folks are using with open source GIS. Ubuntu, some other form of Linux, etc. Thoughts and experiences appreciated.


texts:
I have created a road style with a thinner line overlaying a wider line, with different road classifications on different layers. I've used the advice here to make segmented lines continuous.

This does not work where roads meet:

I wonder if there is a solution to this eg should I take a different approach to rendering? I tried merging all the layers to one then classifying by road type but could not figure out how to then apply the thinner/wider line.

It's OS open data.


texts: Fire department coverage (voronoi equivalent)

texts:
I have a dataset of fire stations within an area. I would like to divide that area up into segments whereby each fire station covers a segment. This should be based on driving time. e.g. if I have 2 fire stations, fire station 1 should attend to all incidents with X segment as it will always be closer than fire station 2.

In the example below, fire station 1 = red dot and 2 = yellow dot. This is fictional, however you will see that station 2 covers a greater area to the south due to the high speed road going west from Hagworthingham.

I have around 50 stations dotted on the map. I'd like to break up their areas accordingly. I can use OSRM or GraphHopper but in reality its the theory behind how I can calculate this that I am trying to understand.


texts: Google earth deviation from arcgis 10 based maps

texts: I really need your help.. I am using arcgis 10 and i mapped some points. with the based maps provided, i tried digitizing the road network. When i exported to kmz file and opened it in google eath. there is a shift. is it that i am not doing something right


texts:
I am using Geopackage stored in a folder using Onedrive.
Everytime I want to modify my file in QGIS (3.12), it is freezing and onedrive is notifying errors in synchronisation for the geopackages files (gpkg-shm + gpkg-wal).
When I use shapefiles, I don't have any problems.

Does anyone faced the same issue or should I modify some parameters in QGIS ?


texts: Converting multiple CSVs into shapefile using Python for loops

texts:
I have hundreds of CSV files with the naming convention similar to:
prcp_025_01_01.csv
prcp_025_05_01.csv
prcp_025_08_01.csv
prcp_050_01_01.csv
prcp_050_05_01.csv
prcp_050_08_01.csv
prcp_025_01_15.csv
prcp_025_05_15.csv
prcp_025_08_15.csv

I was wanting to create a shapefile for each CSV like so:
prcp_025_01_01.shp
prcp_025_05_01.shp
prcp_025_08_01.shp
prcp_050_01_01.shp
prcp_050_05_01.shp
prcp_050_08_01.shp
prcp_025_01_15.shp
prcp_025_05_15.shp
prcp_025_08_15.shp

When running the following code, it will create the first shapefile of
prcp_025_01_01.shp
But then error out with

ERROR 000725: Layer Name or Table View: Dataset prcp_dly_layer already
exists. Failed to execute (MakeXYEventLayer).

Im running 10.5.1
Below is my script. Any suggestions on how to fix it?
import arcpy
import os
arcpy.env.workspace=r"C:\path"
shps=arcpy.ListFeatureClasses()
SR=arcpy.SpatialReference(4269) #WGS 1984

for shp in shps:
    print arcpy.Describe(shp).spatialReference.name

elem = ["025", "050", "100", "400"]
month = ["01", "05", "08", "12"]
day = ["01", "15", "31"]
for element in elem:
    for mon in month:
        for dy in day:
                in_table=r"C:\path\csv_files/prcp_" + element + "_" + mon + "_" + dy + ".csv"
                prcp_dly_layer=arcpy.MakeXYEventLayer_management(in_table, "LON", "LAT", "prcp_dly_layer", SR)

                fc=arcpy.CopyFeatures_management(prcp_dly_layer, r"C:\path\prcp/" + "prcp_" + element + "_" + mon + "_" + dy + "_" + '.shp')

print(arcpy.GetCount_management(prcp_dly_layer))


texts: Data format for checking a land survey deed

texts:
A commercial property is in the closing phase of sale and the only hold up is that the deed description is one that's been lifted from prior deeds dating back to the 1950's and is obviously erroneous.  A 1992 survey does appear to properly depict what the Seller's have to convey and both Buyers and Seller agree that the aforementioned survey be the basis of description in the new deed and title policy if ...it plots out in reasonable fashion in Copan.

But, after spending six hours trying to get a result (closure error, or an additional call necessary to close) and a graph, I have been unsuccessful.
format I've used is:

Map Traverse
1 degrees distance tab tab tab
2 degrees distance tab tab tab
3 degrees distance tab tab tab
1 degrees distance tab tab tab

That produced a bizarre graph, so I tried making the first call

Map Traverse
1 0.0 0.0 tab tab tab

Then I get an error message saying in effect "insufficent tabs" and use cut and paste to provide missing tabs.

My data:

1
2   354.0847    129.65
3    84.4149    117.82
4   352.2410     69.78
5   309.4352     82
6    37.1300     55
7   307.4700     61
1   276.5342     95.62

The original survey plat:


texts: Dump a GeoJSON FeatureCollection from Spatialite query with Python

texts:
Given a Spatialite database, holding a table with a geometry column and some other attribute columns. I want to query that table in a Python script and get a GeoJSON FeatureCollection string in return which stores one GeoJSON Feature per result row, each feature being composed of the geometry from the geometry column and all the other columns as attributes.

Connecting to a SQLite database with Spatialite extensions and doing a query in Python:

import sqlite3

# open database file
conn = sqlite3.connect('MyDatabase.sqlite')

# load spatialite extensions. make sure that mod_spatialite.dll is somewhere in your system path.
conn.enable_load_extension(True)
conn.execute('SELECT load_extension("mod_spatialite.dll")')

# get a cursor
cursor = conn.cursor()

# compose a query
query = """
SELECT
    geometry,
    firstAttribute,
    secondAttribute,
    thirdAttribute,
FROM
    MyTable
;
"""

# execute query
cursor.execute(query)

# get results
results = cursor.fetchall()

However, results is now a list of tuples (one tuple per row) filled with strings (I never understood why this could be useful ever...), i.e.

[(b'\x00\x01\xe6\x10\x00\x00\xde\xe4\xb7\xe8d\xa9\x95\xbf\xf9\xf5Cl\xb0\xbcI@\xde\xe4\xb7\xe8d\xa9\x95\xbf\xf9\xf5Cl\xb0\xbcI@|\x01\x00\x00\x00\xde\xe4\xb7\xe8d\xa9\x95\xbf\xf9\xf5Cl\xb0\xbcI@\xfe', 'value1', 'value2', 'value3'), ...]

The geometry column is given as binary (as stored in the database), however, getting it as GeoJSON representation would be possible by changing the query from SELECT geometry ... FROM ... to SELECT AsGeoJSON(geometry) ... FROM ....

Still, this is miles away from the format I want to achieve in the end.
Any suggestions how this could be done elegantly without having to split the tuple into smallest parts and recomposing them manually to the desired format?


texts: Split feature class based on number of records arcpy

texts: I am looking for a way to split a feature class into multiple feature classes based on the number of records, probably using python. E.g first 100, next 100, etc. I have so far only seen options to split by attributes.


texts: I have a .tif file I would like to use as my base map on the Android version of QGIS. However, when I try and load it as a raster layer, it states that this is not a correct raster source.  Is there a way around this?


texts: How to copy selected features from one sde feature class to another with different attributes using Arcpy?

texts:
I am new to ArcPy scripting, I need help on copying the features to SDE feature class.

My Scenario: We have two Geodatabases running on Oracle 11g ArcSDE 10.2. I have been asked to write the Python script for query and copy the features from SDE connection1: Feature class1 to SDE connection2: Feature class2.

The feature classes are not identical. I need to copy Shape and only two attribute values feature class1 to feature class2 and also populate the remaining attribute values of feature class2 with constant values.

Please guide me in right direction to achieve this.

My Code till now:

# Import system modules
import arcpy
import os
import sys
try:
        ParcelSAMNo = arcpy.GetParameterAsText(0)
        ModeofQuery = arcpy.GetParameterAsText(1)
        TransactionNo= arcpy.GetParameterAsText(2)
        arcpy.env.workspace = r"C:\Users\chenna.kishore\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\LIC.sde\LIC.LandPlan"
        arcpy.env.overwriteOutput = True
        fc = "LIC.LPLN_CADASTREPLOT"
        where_clause = "PIN in (" + ParcelSAMNo + ")"
        s_cursor = arcpy.da.SearchCursor(fc, ("SHAPE@", "PIN", "PDAREA"), where_clause)
        arcpy.env.workspace = r"C:\Users\chenna.kishore\AppData\Roaming\ESRI\Desktop10.2\ArcCatalog\SDE@REGEO.sde\SDE.Parcels"
        outFc="SDE.RPAR_PD_Parcels"
        with arcpy.da.InsertCursor(outFc,["PIN", "PAR_AREA", "SAM_NO", "TRANS_NO", "SHAPE@"]) as rowInserter:
            for s_row in s_cursor:
                    # Insert the new row into the Parcel Layer
                    geom = s_row[0]
                    row_values = [str(s_row[1]), s_row[2], 0, int(TransactionNo), geom]
                    print row_values
                    rowInserter.insertRow(row_values)
                    #rowInserter.insertRow([str(s_row[1]), s_row[2], 0, int(TransactionNo), geom])
        # Clean up the cursor
        del rowInserter
        del s_cursor
        del s_row
except Exception as e:
    print e.message
print "Success Full..!"

I cannot view the data that has been copied. But there is no errors from the script. Any Suggestions?


texts: QGIS - Split a line having attributes with UNIQUE constraint

texts:
I am using QGIS with a PostgreSQL/PostGIS database, where the lines table (named 'segments'), besides the integer serial primary key, has a UUID column having a UNIQUE constraint. Trying to split a line with the SplitFeatures tool of the Advanced Digitizing toolbar fails to commit to the database because the new line that is created after the split, has the same uuid of the original line and therefore violates the database constraint.

In order to simulate the way the primary key is generated using the PostgreSQL sequence generator

nextval('segments_segment_id_seq'::regclass)

, I have tried to manually set a default value in the database using a random uuid generated by:

uuid_generate_v5(uuid_ns_url(), ( (md5(((random())::text || (clock_timestamp())::text)))::cstring))

but this does not work. QGIS continues to duplicate the uuid of the original line to both lines that occur after the split. (Tested that after removing the UNIQUE constraint)

Also tried to include the uuid column into the primary key definition, along with the integer serial. Did not work either.

I assume that QGIS must have a way to identify such columns with constraints when splitting features, because it handles properly the auto-increment of the primary key.

Environment: QGIS 2.8.5-Wien, PostgreSQL 9.4


texts: Using iteration in Model Builder for Clipping features in an ArcGIS File Geodatabase

texts:
I have to Clip the data in a file Geodatabase in which there are about 80 Feature Datasets which includes about 20 feature classes in each one of the Feature Datasets.

I have created the model once which iteration of "Feature Classes", once with "Datasets" and once with " Workspace".

I was able to do some clipping but None of them do it through the whole GDB and just get the first Feature dataset. How can I make it automatically for the whole FGDB?


texts: PyQgis overwrite shapefile and adding layer to QGIS

texts:
I am creating a shapefile in Python by:

driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.CreateDataSource('c:/temp/toke2.shp')
layer = datasource.CreateLayer('layerName',geom_type=ogr.wkbLineString)

If the file toke2.shp doesn't exist, it will be created. If it exists, I can't run this code and overwrite the old file. I need an overwrite function or a delete function. I run QGIS in admin mode.

My other problem is when I want to add the shapefile to QGIS by:

vlayer = QgsVectorLayer("c:/temp/toke2.shp", "toke2", "ogr")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

If I add the file in QGIS by "Layer > add vector layer", it works fine. When I use the code in the python console in QGIS it also works fine.
But if I use this code in my plugin the file is added to QGIS, but it is empty. So for some reason it doesn't add the data points in the file, it is just adding an empty shapefile.

In my plugin I import this:

 from osgeo import ogr
 from osgeo import osr
 import math
 from PyQt4.QtCore import *
 from PyQt4.QtGui import *
 from qgis.core import *
 import qgis.utils
 import resources_rc

Does anyone have a solution to these problems?


texts:
I have a list of coördinates that are in X/Y form, and I need them in lat/lon.
I tried googling around a bit, and found this method:
Converting projected coordinates to lat/lon using Python?

This does seem to work, but I found out a few minutes ago that there are numerous coordinate reference systems, and I have no clue which one my data is in.

I tried a few of these: https://epsg.io/?q=Belgium (because the file I got is from a Belgian company), but the results I got didn't seem right.

Is there any way to figure out what coordinate system was used (except for trying out all 10000 possibilities), with a known combination?

For example, I know X="0511245565" Y="0042422978" should approximately equal to 51.206728, 3.253426 in wsg84.

My appologies if this has been asked before, or if this is a dumb question. I have no background in GIS, or anything related.


texts: Can gee geometries be created using server-side variables?

texts:
Say I am calculating some numbers in some way and I end up needing to create some geometries based on those numbers. The simplest example would be:
var bboxList = ee.List([1,2,8,9]) // assume this has been calculated in gee
// the below line should, in my opinion, work, but it doesn't,
// because ee.Geometry.BBox can only handle client-side objects
var myBbox = ee.Geometry.BBox(bboxList.get(0),bboxList.get(1),bboxList.get(2),bboxList.get(3))
// what finally works is getting the server-side object to the client-side,
// but this is highly discouraged
var myBbox = ee.Geometry.BBox(bboxList.get(0).getInfo(),bboxList.get(1).getInfo(),bboxList.get(2).getInfo(),bboxList.get(3).getInfo())

Is there a way to create geometries, without passing all their coordinates by the client?
It is obviously not an issue when there's only 4 coordinates, but one can easily imagine this becoming an issue with thousands, or millions of coordinates sent back and forth.


texts:
I'm having problems with Molusce pluging. The rasters have same cell size, coordinate sytem and extension. My problem is that I get the transition Matrix but the information,  according my experience is not correct. The highest values should be in the columns with the same land cover that in the files.

This picture shows that the urban land cover has the highest values in Suelo Desnudo and I'm almost sure that isn't correct. The goal is to know in which land cover has made the transition every land cover.

Why is this happening?


texts: Linear fit in GEE, why is the ouput plotted as an RGB image and what do the max and min values signify?

texts:
I am currently reading up on the linearFit() and linearRegression() in GEE. I was looking at the developer guide for these functions, which I think I understood correctly, except for the way in which the outputs are added to the map. Once the linear regression or fit is performed the outputs are plotted as an RGB image assigning R and B to 'scale', and G to 'offset', like this:
Map.addLayer(linearFit,
  {min: 0, max: [-0.9, 8e-5, 1], bands: ['scale', 'offset', 'scale']}, 'fit');

According to the guide the resulting map shows areas of increasing trend in blue, decreasing trend in red and no trend in green, which seems simple enough but I don't understand how the max and min values were chosen, or the point of using and RGB composite is, versus, for example, just plotting 'scale' with any given palette.


texts: SDO_GEOM.RELATE with EQUAL mask returns false for same WKT