texts: excluding NaN records from FOR loop in ArcPy?

texts:
I am trying to create an array of Z and M values for vertices from a polyline (Z and M enabled) feature class.  Some of the vertices do no have M values/geometry and I wish to exclude them from the array.

    import arcpy
    from arcpy import env
    ##import math   #another way I was testing....commented-out
    import numpy as np

    print "Modules imported."

    # Set environment settings
    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = r"E:\Projects\HaulRd\SpatialData"
    arcpy.env.scratchWorkspace = r"E:\Projects\HaulRd\SpatialData\Scratch"
    arcpy.env.outputZFlag = "Same As Input"
    arcpy.env.outputMFlag = "Same As Input"

    print "Workspaces and environments set."

    # Set local variables
    infc = "HaulRd_XS_linesZM_test.shp"
    ##x=float('nan') # another way I was testing.....commented out

    # Identify the geometry field
    #
    desc = arcpy.Describe(infc)
    shapefieldname = desc.ShapeFieldName

    # Create search cursor
    #
    rows = arcpy.SearchCursor(infc)

    # Enter for loop for each feature/row
    #
    for row in rows:
        # Create the geometry object
        #
        feat = row.getValue(shapefieldname)

        # Print the current feature's ID
        #
        print "Feature %i:" % row.getValue(desc.OIDFieldName)
        partnum = 0

        # Step through each part of the feature
        #
        for part in feat:
            # Print the part number
            #
            print "Part %i:" % partnum

            # Step through each vertex in the feature
            #
            for pnt in feat.getPart(partnum):

                #exclude records with 'nan' (i.e. null)
                if pnt[~np.isnan(pnt).any(1)]

                    # Print z,m coordinates of current point
                    #
                    print pnt.Z, pnt.M

                else:
                    # If pnt is None, this represents an interior ring
                    #
                    print "Interior Ring:"
            partnum += 1

I have made a number of attempts based on content I found in:

https://stackoverflow.com/questions/11620914/removing-nan-values-from-an-array
https://stackoverflow.com/questions/2695503/removing-pairs-of-elements-from-numpy-arrays-that-are-nan-or-another-value-in
https://stackoverflow.com/questions/12986996/conditional-for-in-python

But I can't seem to figure-out how the syntax (could also be a logical problem) fits into my code.

Eventually, my temporary array will be be used to feed creation of vertical line graphs (see related post:

Using Z and M feature geometry to create line graph in ArcPy/Python?)

I am using Python v2.6.5.


texts: Can I access recent formula in qgis field calculator?

texts: I'm working with QGIS field calculator and I sometimes need to adjust the formula for a field a little bit. My problem is, that whenever I adjust the formula, I need to type the formula as a whole again. When using the field calculator in ArcGIS the formula is not lost after executing it, is there a way that's possible in QGIS? It costs me so much unnecessary time and effort...


texts: How to get WMS layer title instead of layer name from Geoserver?

texts:
I specified the raster file with "name" and "title" in Geoserver.

For example, I have these in Geoserver:

name,  Title
AT:usa, "This is a map of USA"

Can I return the title ("This is a map of USA" ) as the layer name in Openlayers layer switcher when I call the WMS service for this layer by name? By default, it shows "usa" in the layer switcher.


texts: Qualifying field names in field calculator

texts:
I am trying to set up a quick model that:

Takes a point feature layer of development projets (a table in an access database with xy coordinates),
Spatially joins polygon neighborhood designations to the feature layer,
Join the output of 2 to the source feature layer (sitting as a table in access)
Updates the neighborhood field based on the neighborhood field from the attributes from the spatial join

The point is that I really don't want to create a new file every time I do this, but maintain and add to the existing table of development projects.

My question is how to properly qualify field names in the field calculator so I can update values in the source table (1) based on the values of the spatial join featureclass with the neighborhood designation (2). I join the output of (2) to (1) so the attributes should all be in the same name space, I imagine.

I currently use this as my update expression (VB parser):

[PipeRealCurrent.PlanningDistrictsCombo]=[PipeRealCurrent_x_PDist.PlanningDistrictsCombo]

A field name was not found or there were unbalanced quotation marks.
ERROR 999999: Error executing function.
Item not found in this collection.
Failed to execute (Calculate Field).


texts:
Recently I started learning Openlayers 2x and I've been struggling for a while now with a GeoJSON and it's projections.

So, here is the deal;

I successfully added wms layers from ( http://geoportal.dgu.hr/podaci-i-servisi/svi-servisi-i-aplikacije/ ) with a projection HTRS96/TM (EPSG:3765) and then when I try to add vector layer it doesn't render it. Vector layer was created using QGIS software with a projection code:

"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3765" }

Here is the  OpenLayers code:

var map;

function init() {
    map = new OpenLayers.Map('map_element', {
        projection: 'EPSG:3765',
        maxExtent: new OpenLayers.Bounds(426906.034398, 5101783.64664,
        468058.881494, 5123976.6968),
        numZoomLevels:15,
        maxResolution: 'auto',
        units: 'm'
    });
    var wms = new OpenLayers.Layer.WMS(
        'DOF 1:5000',
        'http://geoportal.dgu.hr/wms?layers=DOF',
        {}
    );

    var wms2 = new OpenLayers.Layer.WMS(
        'TK 1:25 000',
        'http://geoportal.dgu.hr/wms?layers=TK25',
        {}
    );

        var wms3 = new OpenLayers.Layer.WMS(
        'HOK',
        'http://geoportal.dgu.hr/wms?layers=HOK',
        {}
    );

    map.addControl(new OpenLayers.Control.LayerSwitcher({}));
    map.addControl(new OpenLayers.Control.MousePosition());

    var KZZ = new OpenLayers.Layer.Vector("KZZ", {
    protocol: new OpenLayers.Protocol.HTTP({
        url: "KZZ_JLS.json",
        format: new OpenLayers.Format.GeoJSON()
    }),
    strategies: [
        new OpenLayers.Strategy.Fixed()
    ]
});

    map.addLayers([wms,wms2,wms3,KZZ]);
    map.zoomToMaxExtent();
}

And here is GeoJSON layer (fragment of it, whole file can be downloaded from this link https://dl.dropboxusercontent.com/u/126827789/KZZ_JLS.json):

{
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::3765" } },
"features": [
{ "type": "Feature", "properties": { "POVRINA": 27.027842695, "NAZIV": "Stubièke         Toplice" }, "geometry": { "type": "Polygon", "coordinates": [ [ [ 457191.01612696185, 5094147.6274271049 ], [ 457199.98992743361, 5094130.2420611409 ], [ 457231.34151837695, 5094112.0906051407 ], [ 457283.71295271069, 5094094.2316290718 ], [ 457310.09355065826,....

How can I render GeoJSON and add it to this file?

This problematic openlayers map can be viewed here htt()p://interaktiva.pregrada.hr/nc/--> to low reputation to add 2nd link on my post, pls remove ()


texts: ESRI Javascript API: Return to the Map Extent you were at when you use browser back button

texts:
Can anyone give me some pointers (code would be helpful) on how to go back to the extent a user was at after they have navigate away and clicked the back button?

So I'm viewing a map at X extent I click on some link and page changes, on clicking the back browser button I want to be able to get back to the Map Extent I was just before the page changed.

Any help much appreciated.


texts: PostGIS geocoding with local street centerline data

texts:
I have the PostGIS TIGER geocoder setup and working to geocode addresses in a single city. I am trying to use local street centerline data for geocoding instead of the TIGER data because it is more complete.

The linked PostGIS documentation briefly mentions that the geocoding function could be altered and applied to other data sets. What approach would I need to take to alter the PostGIS geocoding function to use a different data set?


texts: Is there a way in the Google Maps API to perform a union?

texts: I'm looking for a way to take the overlap area of multiple polygons and work with that result as a single entity, ideally all in the API.


texts: Workpath for Tiled GeoTIFF Warp and Worldfile

texts:
I am attempting to identify a workpath (or script) for warping TIFFs (and TFWs) and saving them out.
Have tried diff methods "batch mode" and pointed it to the relevant folder etc. But I just get an error message. See pic attached.

Is there any settings that are incorrect? Or a script I could run that does the whole routine?

JPEG with JGW is the ideal output. But not essential, where warping is.


texts:
I have a PostgreSQL database with multiple schemas. I am trying to save data using Arcpy to my database but I cannot figure out how to define the schema I want it in.

The postgres.sde database has 3 schemas; how do I differentiate where the table will be saved?

Here is the test performed:

input_fc = X.shp

output_location =

(TRY1 -
'C:\\Users\\david\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\postgres.sde'),

(TRY2 - 'C:\\Users\\david\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\postgres.sde'),

(TRY3 - 'C:\\Users\\david\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\postgres.sde\\projects\\david')

output_table_name = (TRY1 -'a3_table'), (TRY2 - 'projects.david.a3table'), (TRY3 - 'a3_table)

RESULT:

TRY1- No error Just a seemingly random location of table output

TRY2- ERROR 000706: Output name projects.david.a3_table cannot contain dots (.)

TRY3- ERROR 000732: Output Location: Dataset C:\Users\david\AppData\Roaming\ESRI\Desktop10.3\ArcCatalog\postgres.sde\projects\david does not exist or is not supported

arcpy.TableToTable_conversion(input_fc, 'C:\\Users\\david\\AppData\\Roaming\\ESRI\\Desktop10.3\\ArcCatalog\\postgres.sde', 'a3_table')


texts: I am curious if anyone has come up or has a good graphic style they can share that is symbology for NAICS codes. I have a large dataset that I want to symbpolize based on the 2 digit NAICS industry code, but can't any sort of good or consistent stanard.Any good sources?


texts: Questions on smooth lines on soil polygon?

texts:
In a couple of weeks I am going to do QA/QC on soil polygons and I'm wonder if there is a way for me to run a script or a tool box that allow me to find the jagged lines ?

Alternative, my plan would be to set up using a quad layer for each quad and then review the lines that need to be smooth. I would create a new layer to mark an error that need to be correct smooth line.  I feel that it would take about a week or so to complete the work. The layer itself is a big file.

In the past, I have had convert from polygon to polyline but it took a lot of time and work to clean up. If you do that, there would be probably a lot of gaps in between lines ane would create gaps and produce topology errors.

Does anyone have a suggest a better way around it ? I am using ArcGIS 9.3.1 SP3.

Thanks !


texts: How to display only one placemark on Google Maps via JavaScript?

texts:
I have a KML file which contains multiple placemarks. I'm displaying this file on Google Maps following this description: https://developers.google.com/maps/documentation/javascript/kml

Is it possible, using the JavaScript API, to select only one placemark to be displayed?


texts: Reading GeoJSON-String from CSV as Geometry in QGIS

texts:
I have a CSV which contains the geometry information as a GeoJSON string, e.g.
{ "type": "MultiLineString",
    "coordinates": [
        [[10, 10], [20, 20], [10, 40]],
        [[40, 40], [30, 30], [40, 20], [30, 10]]
    ]
}

which is just a text field.
I would like to display this geometry in QGIS 3. I could go via PostGIS and use ST_GeomFromGeoJSON, but would prefer a direct way, without adding my data to a PostGIS database first.
Unfortunately, virtual layers do not recognize this function and as expressions there only seem to be geom_from_wkt() and geom_from_gml().
Is there a possibilty to read this GeoJSON-String as Geometry directly from a CSV column?


texts: Geopandas shows points outside of buffer geometry even though inside

texts:
I have some ship trajectories and I want to know which point of the trajectory is inside the port(=buffer) and which points are outside. The following picture shows my Geopandas Dataframe with the geometry.

The problem is that some points that are inside, are shown as outside. The problem can be easily seen if I use a huge buffer.

After that I am loading the port's point and buffing it, like below.

And then I am finding which points are inside and which outside, with the code below.

Now if I use a huge buffer, the are still rows in the points_outside_geometry dataframe. These rows shouldn't be tracked as outside points and this can be verified by the following plot.

How is that possible?

Code:

polys = gpd.GeoDataFrame.from_file('final_port.shp').buffer(100000)
sindex = gdf.sindex
points_within_geometry = pd.DataFrame()

for poly in port_poly:

    possible_matches_index = list(sindex.intersection(poly.bounds))
    possible_matches = gdf.iloc[possible_matches_index]
    precise_matches = possible_matches[possible_matches.intersects(poly)]
    points_within_geometry = points_within_geometry.append(precise_matches)

points_within_geometry = points_within_geometry.drop_duplicates(subset=['MMSI', 'TIMESTAMP'])
points_outside_geometry = gdf[~gdf.isin(points_within_geometry)].dropna(how='all')


texts: Error: QGIS standalone script in crontab

texts:
I am trying to automate QGIS (2.18.4) standalone scripts in ubuntu using crontab.
My scripts execute well when running from command, but fail if I put in crontab.

from qgis.core import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.gui import *
app = QApplication([])
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

In crontab, it fails at app = QApplication([]). Crontab gives error as

  : cannot connect to X server

How to make my QGIS standalone script run in crontab?


texts: Calculating distance between two points on map using ITM coordinates?

texts:
How to calculate the distance between two point on the map using ITM coords?

i have a map with items on it all with ITM coords and i need to return all items in X meters of radius from my location.

the problem is how to convert the distance between the two points on the map to meters or km.


texts: Live Tracking with Garmin inReach explorer+

texts:
I'm trying to connect my Garmin Inreach Explorer+ gps to QGIS (running on a MS Surface 5 tablet) via bluetooth, but have been unsuccessful so far.  The goal is to have to a live tracking tablet while working in the field, passively collecting GPS coords while mapping with QGIS on the tablet.  The GPS pairs and connects to the tablet just fine, however when i try to connect with the live tracking module on qgis, it doesn't work.  Not too sure what else to do...

I have also tried downloading GPSGate, but this also didn't work.  Is it possible to connect a Garmin gps into QGIS for live tracking in the field?


texts: Comparing real orthophotos against common available tiles on the internet (Google, Bing, ..)

texts:
All available background tiles I know about on the internet (OpenLayers, Bing, Google) are tilted. See a random example below from OpenLayers (you can see the east side of the buildings):

In essence, when you move the mouse from left to right along the same building, the values of the XY coordinates change, how is that?, should not the coordinates remain the same while you are on the same 2-D position?. Obviously OpenLayers does not correct that, neither the others, am I right?

Is there any way of taking this effect into account? Basically what I want to do is overlapping a true vertical digital terrain model (DEM) with aerial tiles on top. Should I worry about this distorsion?


texts: How to change label of rotation button in OpenLayers?

texts:
The documentation claims that "instead of text also an element (e.g. span) can be used":
https://openlayers.org/en/latest/apidoc/module-ol_control_Rotate.html
I tried
this.olMap = new Map({
  // replace arrow with compass icon on rotation button
  controls: defaultControls({
    rotateOptions: {
      autoHide: false,
      className: 'my-compass',
      label: '<svg class="ol-compass compass-svg block w-5 h-5 cursor-pointer select-none z-120" role="button" aria-label="Compass"><use xlink:href="#compass" /></svg>',
      tipLabel: 'Reset GeoMover Rotation'
    }
  }),
[...]

However, that does not work. Neither does using the class name "compass-svg" or simply the  tag.
What is the proper use of label with a svg element?
UPDATE:
In addition to the accepted solution, I'd like to add an example for incorporating a "use" fragment with svg, which needs special care when generated with Javascript (see https://florianbrinkmann.com/en/svg-use-element-javascript-4513/):
<!-- HTML -->
<!-- A svg element that only "collects" symbols: -->
<svg xmlns="http://www.w3.org/2000/svg" style="display: none;">
  <def>
    <!-- This fragment shall be used for the rotate button -->
    <symbol id="compass" xmlns="http://www.w3.org/2000/svg" fill="none" viewBox="0 0 24 24"><g clip-path="url(#clip0)"><path fill="currentColor" d="M10.98 2.19L7.5 11.48a2 2 0 000 1.4l3.47 9.29a1.24 1.24 0 002.31 0l3.47-9.3a2 2 0 000-1.39l-3.47-9.3a1.24 1.24 0 00-2.31 0zm3.15 9.99a2 2 0 11-3.99 0 2 2 0 014 0z"/></g><defs><clipPath id="clip0"><path fill="#fff" d="M0 0h24v24H0z"/></clipPath></defs></symbol>
  </def>
</svg>

// JavaScript

// function that returns the svg element
compassIcon () {
  const svgElem = document.createElementNS('http://www.w3.org/2000/svg', 'svg')
  const useElem = document.createElementNS('http://www.w3.org/2000/svg', 'use')

  svgElem.setAttribute('role', 'button')
  svgElem.setAttribute('aria-label', 'compass')
  svgElem.classList.add('ol-compass', 'compass-svg')

  // `#compass` is the id of the svg symbol
  useElem.setAttributeNS('http://www.w3.org/1999/xlink', 'xlink:href', '#compass')
  svgElem.appendChild(useElem)

  return svgElem
}

// Usage with OpenLayers in Vue
export default {
   data: () => ({
      olMap: null,
   }),
   mounted() {
      this.olMap = new Map({
          // add rotate functionality
          interactions: defaultInteractions().extend([new DragRotateAndZoom()]),
          // change rotate button
          controls: defaultControls({
              rotateOptions: {
                  autoHide: false,
                  label: this.compassIcon()
              }
          }),
      }
   }
}


texts: Is it possible to export GeoTIF current view to GeoTIFF with QGis?

texts: ArcGis offers this, but I didn't found such a function in Qgis. Maybe I have simply not found the command.... Does someone know the command or is there an extenstion which offers such a functionality?


texts: Can We build a search function using GeoExt?

texts:
Or do I need another technology?

I need to build a web gis server and I need map display, layers switch, search function and attribute display feature.

I am currently using GeoServer, OpenLayers and GeoExt.


texts: Simplifying multiple lines to create central axis?

texts:
I have a layer with railway lines, downloaded from OpenStreetMap. Every single track is represented separately, so at trainstations, you have several parallel lines, as well as sidelines that lead to depots. Sometimes, lines are broken (it's OSM!). I want to reduce all these lines to one single, simple line.

An idea is to create buffers around the railway lines and to dissolve them to a single polygon (and maybe apply a negative buffer to reduce its width and generalize the geometry)?

If I could create a medial axis from this polygon,, it would represent quite well the simplified railway line I want to create.

However, how to create such a medial axis in QGIS?

There is a solution in PostGIS: https://postgis.net/docs/ST_ApproximateMedialAxis.html

But unfortunately, I have no experience with postgis, so I'm looking for a solution inside QGIS. A possible solution seems to be the HCMGIS-plugin that has a function "Create Skeleton/Medial Axis". The result represents quite well the railway-track (see screenshot, red line), but with one major problem: it creates a skeleton with small side-lines connecting to the vertices of the simplified buffer.

Second problem: algorith takes quite long for a larger selection - I was able to test it only with a small selection as the calculation for the whole railway-network did not stop even after half an hour. I'm not sure if it works well for the whole network (see screenshot 3 with the red part representing the selection used for testing).

Any idea for alternative solutions?

@Cyril's suggestion of using the midline of the skeltonize method would be perfect for me, but unfortunately, the resulting layer consists of only one object - so the short lines are not separated from the long one.

Maybe there is also a completely different method for solving my problem?

Update: see PostGIS-syntax: export query-result to QGIS / Shapefile


texts: Default Opengeo suite print module- scales given in selector differ when using the zoom slider

texts:
Using a default installation of the 3.0.2 suite, opening the Composer view(with the MapQuest OSM tiles as background) and clicking print (to pdf) to activate the print preview there appears to be a couple of issues with the scale given in the scale selector box.

One- if one slides the zoom slider so that the scale selector box gives 1:512 000, and then opens the scales list in the box and chooses 1:512 000 the scale visibly changes in the print preview map (it goes to the next smaller scale).

Two- the scales in the scales list in the scale selector box do not entirely match the list of scales given in config.yaml: the smallest (1:524288000) is not listed.

Three - there are 17 zoom levels in the default Composer map, 18 in the print preview's zoom slider, and 19 in the scale selector box.

The scales given by the zoom slider in the default Composer map roughly match the scales given in Config.yamls, i.e 1:68k something becomes 1:64k etc, but this should not be reason for this default behaviour?

I have tried changing the number of scales config.yaml as well as changed the maxZoomLevels for the map, but besides zoom slider in the print preview following the number of zoom levels (scales) of maxZoomLevels, even if the main map's zoom slider does not, I cannot get things to precisely line up.

A search here and abroad did not reveal this issue described, but if parts of this (or all) is due to a lack of understanding of certain fundamentals on my side, then pray tell where I can go read :)


texts: Looking for simple software to compare two Landsat Images

texts: At uni I used ERDAS imagine to view two rasters at the same time and it had various viewers for comparing them.  I'm looking for something that has a slider so that I can see the differences between two Landsat images (same area different dates) or if not a slider some good functionality for visually comparing them.  ArcMap is too slow at switching between the two.  Is there anything free that I could use for this?


texts: I'm using the OSM map and shp. file to create a map with a rive but I can't get a grid up in the layout manager and have tried everything.. The projection is set by default to IRENET95 / Irish Transverse Mercator but the OSM map is WGS 84 / Pseudo-Mercator.. could this be the issue? Thanks!


texts: Select tracks or part of tracks that are within a polygon?

texts:
I am studying a monkey species, I have GPS observation points which I used to make a convex polygon (Hull). I also have all of the survey trails imported and displayed in QGIS (See screen shot for clearer understanding). As I have over 2700hr of survey tracks I want to select only the tracks or parts of the tracks that are within this polygon. The end goal is to get a more accurate estimate for survey effort within the core range (convex polygon). (using QGIS 2.18.2)

I have tried vector->research tools-> select by location which only found tracks that were wholly within the convex hull shapefile.

I then tried vector->geoprocessing tools -> clip which wasn't able to run giving the following error:

2017-05-16T15:45:56 2   Uncaught error while executing algorithm
            Traceback (most recent call last):
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/core/GeoAlgorithm.py", line 203, in execute
                self.processAlgorithm(progress)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/algs/qgis/Clip.py", line 71, in processAlgorithm
                source_layer.crs())
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/core/outputs.py", line 326, in getVectorWriter
                crs, options)
              File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing/tools/vector.py", line 610, in __init__
                uri = GEOM_TYPE_MAP[geometryType] + "?uuid=" + unicode(uuid.uuid4())
            KeyError: 1005

FINAL EDIT

Once converting my track file from multiparts to single parts (Vector->geometry tools-> multiparts to single parts). This was difficult as a newbie as most people refer to it as MULTILINESTRINGZ and LINESTRING on stackexchange. I was then able to clip the track shapefile to the convex hull polygon.


texts: I have 2 layers in OpenLayers: first is parcel layer (polygons) and second is railway layer (lines). They are already drawn and all I need to do is when I click on some line, to get intersection of that line and all polygons that are intersected with selected line and get feature info of polygons. How can I do that?


texts: The best way to interpolate for regular timed data

texts: I have 20 stations, in the table of which I have the coordinates and the value of the other 3 factors for 2020 with regular intervals of one day. Now I want to create a  raster for each of these 3 factors and for each day, which is the result of interpolation. Is there a way to create this large number (365*3) of rasters faster?


texts:
I realize this is a bit of an open question but are there any known open JavaScript libraries which facilitate printing an OpenLayers map?  I'm looking to produce the content on the client and possibly upload it to convert to PDF or simply set the innerHTML on some DOM element.

I started a proof-of-concept but was looking for something which honors various tile sizes and vector markup (IE 8) and all the other complexities that are sure to arise.


texts: Geocoding building with different entrances?

texts:
I am working on a routing map for a university campus. the routes are modeled and network analyst works fine.

The problem: Many of these buildings have several entrances.

Is there a way to geocode all of them and have network analyst understand which entrance is most suitable or do I have to geocode them with different names?

I ran a test and it seems, the address locator always picks the leftmost entrance with the similar name.


texts: pyqgis: add attributes of points to attributtable of lines

texts:
I'm using QGIS 3.10 with Python 3.7.

I have two layers: "line_layer" and "vertices_layer". I want to add only the names of my vertices (attributtable: "node_name") to the attributtable of my "line_layer" (like start and end node --> attributtable: "node_start" and "node_end"). I have a net so for one line a point A is the start, while for another line this point A is the end. Thats why I can't divide my vertices into start-/endpoints. I want to do this in the console using python (I'm a beginner). My (definetly wrong) code is:

for f in line_layer.getFeatures():
    geom = f.geometry().asMultiPolyline()
    start_point = QgsPoint(geom[0])    #first point
    end_point = QgsPoint(geom[-1])     #last point
    f.['node_start'] = QgsGeometry.nearestPoint(start_point)
    f.['node_end'] = QgsGeometry.nearestPoint(end_point)

line_layer.commitChanges()

The code doesn't work and I don't know how to get only th name of the nearest point in the attributtable of my line_layer.

Just for understanding the situation:


texts: Producing plant height with digitial surface model

texts:
I am trying to produce plant height with a digital surface model "TIF" file that was obtained with a drone. I have clipped out the section with plants and then used "Gap Fill" to produce a digital terrain model. I then used the "Minus" function between the surface model and terrain model.

Does anyone have any other ideas on how to do this that are not so time consuming. The results I got were not very accurate either compared to the in field measurements.

The picture below is of sunflower plants that I would like to get the height of.  I would need the height of an entire section.  Not a single plant.


texts: Is it possible to change the default "symbol template" for layers added to QGIS?

texts:
Can QGIS be configured so that the fill symbols of layers added to the map automatically have a transparent border?  I'm calling it a "symbol template" because I don't know what that's called in QGIS-world.

I would like to set the border to transparent, but keep the random fill color (or however QGIS is choosing the next color) for all layers added to the map.


texts: How to integrate a password-protected WMS layer into an OpenLayers GIS website without entering a separate password?

texts:
I'm building a small intranet site to display geodata easily. Using QGIS (2.18) and QGIS2WEB (2.34.0) I prepared geodata and created it as OpenLayers output as index.html.

A WMS layer, which is protected by username and password, is also part of the output. At the moment, however, an input mask opens, which asks for the username and password. After input, everything works correctly. But I don't want the users of the intranet site to have to enter a username or password.

Unfortunately, the following the proposal, circulating here in the forum, does not work. The input mask still appears.

http://user:password@geoinfo.xxx.de/wms10/72_sd_17_ds_gras.exe

There is no need to hide the data. They can be in the source text.

Is there another way to insert username and password 'automatically'?


texts: Which Image format is recommended/ best for geoserver raster wms services?

texts: I have very high resolution image of vast area(true color RGB bands). I am planning to to reduce the image size by reducing it to 16 bits (or less). which image format is best for performance point of view. finally to narrow down more my question shall i use mbtiles or tiff or mrsid or any other format/ method you like to suggest ?


texts: Change color of feature by ID (GeoTools)

texts:
I am trying to change the color of a specific object, by ID (or another feature). Currently I am able to extract data about an shapefile and display that data. However, I want to be able to change the color of a specific feature (all features that match X ID) based on the ID of the feature (assuming it has one).

Code:

private void changeLineColor() {

    ArrayList tempArray = new ArrayList();

    SimpleFeatureTypeBuilder typeBuilder = new SimpleFeatureTypeBuilder();
    typeBuilder.setName("Features");
    typeBuilder.add("geom-attr-name", MultiLineString.class);
    typeBuilder.add("id", Integer.class);
    typeBuilder.add("id", Integer.class);

    SimpleFeatureType type = typeBuilder.buildFeatureType();
    SimpleFeatureBuilder builder = new SimpleFeatureBuilder(type);
    DefaultFeatureCollection features = new DefaultFeatureCollection();

    try (FeatureIterator<SimpleFeature> iterator = collection.features()) {
        while (iterator.hasNext()) {

            SimpleFeature feature = iterator.next();

            tempArray.add(feature); // all features

        }
    } catch (Exception ex) {

        Logger.getLogger(ShapefileExporter.class.getName()).log(Level.SEVERE, null, ex);

    }finally {

        for (int j = 0; j < tempArray.size(); j++) {
                System.out.println("Attribute: " + tempArray.get(j) );

        }
    }
}

It would look something like this, but only for one MultiPolygon, or MultilineString

I would most likely want to make a dialog like this, but the color would only be for matching codes


texts: Point size based on attribute value in QGIS

texts:
First time user with QGIS. Basically using the MMQGIS plugin I've geocoded a number of points to a map of Ireland. These points reference a club which received a grant. However, I want to be able to distinguish which club received a larger grant.

Is it possible to do this, e.g. various points weighted by grant? So a club which received larger grants have larger points?


texts: I've found yaml files describing GeoServer configuration REST API on docs.geoserver.org website. But where are they on my localhost GeoServer?

texts:
I've found yaml files to generate with Swagger the client that will conveniently call the REST API provided by GeoServer on GeoServer documentation website.
For example at https://docs.geoserver.org/stable/en/api/#1.0.0/datastores.yaml for datastores management :

But I'm expecting to find these yaml somewhere also on my GeoServer installation on localhost (I'm using Geoserver 2.18). They should be exposed at an URL.
But which one? I'm wasn't able to find it until now.

When you go to the yaml documentation, the page that is serving you is using a swagger-ui interface from an internal localhost:8080 GeoServer, and you have even the Swagger "Try it out" buttons.

It means that Swagger knows the structure of the REST services in that localhost:8080 Geoserver.
Aside from this page, which is a perfect and standard one, another one shall exist (that is coming and generated by Swagger too) and that returns these descriptions but in json format, making it compatible for another Swagger tool to generate Client API.


texts: DBF file size is increased after exporting opened shapefile into new shapefile

texts:
I have a shapefile (polygon features) which has seven attributes (id, area, etc.,) in the DBF file. After reading the shapefile I am exporting to a new shapefile in which I am adding two attributes (id ,area). But exported DBF file is bigger size than earlier one.

This is my simple feature type schema

    builder.setCRS(DefaultGeographicCRS.WGS84);
    builder.setName(plotIdentifier);

    builder.add(sPolygon, Polygon.class);
    builder.add(plotIdentifier, String.class);
    builder.add(aRea,Float.class);

    final SimpleFeatureType featureType = builder.buildFeatureType();
    return new SimpleFeatureBuilder(featureType);

Please help me on this

Code creating polygon:

 public boolean addPolygonFeature(List<SimpleFeature> featureList, Point3D[] thePoints,
                                 String polyId, float area)
                                 throws Exception
{
    Polygon aPoly = createPolygon(thePoints);

    if (aPoly == null)
    {
        StringBuffer buf = new StringBuffer();
        buf.append(WError1);
        buf.append(polyId);
        buf.append(WError2);
        String error = buf.toString();

        int opt = GUIBase.showConfirmDialog(error, "Confirm", GUI_YN);
        return (opt == GUIBase.GUI_YES);
    }

    SimpleFeature feature = featureBuilder.buildFeature(null);
    feature.setAttribute(plotIdentifier, polyId.trim());
    feature.setAttribute(sPolygon, aPoly);
    feature.setAttribute(aRea, area);
    featureList.add(feature);
    return true;
}

public Polygon createPolygon(Point3D[] thePoints)
{
    int count = (thePoints == null) ? 0 : thePoints.length;
    if (count < 1) return null;

    // Store Coordinates of polygon as a LINESTRING.
    StringBuffer buffer = new StringBuffer();
    buffer.append("POLYGON((");

    for (int i=0; i<count; i++)
    {
        buffer.append(thePoints[i].x);
        buffer.append(SPACE);
        buffer.append(thePoints[i].y);
        buffer.append(COMMA);
    }

    String lineString = buffer.toString();
    lineString = lineString.substring(0, lineString.length()-1) + "))";
    if (DEBUG) System.out.println(lineString);

    try
    {
        // Creating the Geometry (Polygon)
        return (Polygon) wktReader.read(lineString);
    }
    catch (Exception excep) { }

    return null;
}

public void writeShape(File theFile, List<SimpleFeature> featureList) throws Exception
{
    ShapefileDataStoreFactory dSFactory = new ShapefileDataStoreFactory();

    Map params = new HashMap();
    params.put("url", theFile.toURI().toURL());
    params.put("create spatial index", Boolean.TRUE);

    ShapefileDataStore newDataStore =
                 (ShapefileDataStore) dSFactory.createNewDataStore(params);

    newDataStore.createSchema(featureBuilder.getFeatureType());
    newDataStore.forceSchemaCRS(DefaultGeographicCRS.WGS84);

    String typeName = newDataStore.getTypeNames()[0];
    SimpleFeatureSource featureSource =
                                   newDataStore.getFeatureSource(typeName);
    if (featureSource instanceof SimpleFeatureStore)
    {
        SimpleFeatureStore featureStore = (SimpleFeatureStore)featureSource;
        SimpleFeatureCollection collection =
                  new ListFeatureCollection(featureBuilder.getFeatureType(),
                                            featureList);
        Transaction transaction = new DefaultTransaction("create");
        featureStore.setTransaction(transaction);
        try
        {
            featureStore.addFeatures(collection);
            transaction.commit();
        }
        catch (Exception exp)
        {
            exp.printStackTrace();
            transaction.rollback();
        }
        finally
        {
            transaction.close();
        }
    }
}

And shapefile link as given below:

https://docs.google.com/file/d/0B5Vh2tWJlQW_RURQaFBwaDJ6ckk/edit?usp=sharing

Read shapefile: 1.shp
Written shapefile: 2.shp


texts: Passing data into leaflet js via a php script

texts:
I'm trying to display route data on a map using leaflet.  I'm pulling the data from a MYSQL database with a php script. which is working.  The data is stored in an array called vechiclelatlong.  I use <?php include('../php/vTrackerv3.php'); ?> under the Javascript tag where is have my Leaflet script.  I am getting the following error:

  Uncaught ReferenceError: vehiclelatlong is not defined

See scripts below:

PHP script:

<php
function vData(){
include("dbconnecti.php");

$strtDate = strtotime('20150404');
$endDate = strtotime('20150405');

$query = "SELECT ED.latitude, ED.longitude, ED.timestamp, ED.creationTime
FROM EventData ED
WHERE ED.latitude <> 0 AND ED.longitude <> 0 AND (ED.timestamp >= $strtDate)  AND (ED.timestamp <= $endDate)";
$results = $dbconnect->query($query);

if (!$results)  {
    echo mysql_error();
    die('There was an error excuting query statment: ' ); //. mysqli_error());
}

$num_of_rows = $results -> num_rows;
$data = array();

echo "var vehiclelatlong = [";

for ($x = 0; $x < $num_of_rows; $x++)   {
    $data[] = $results -> fetch_assoc();

    echo "[".$data[$x]['latitude'].",".$data[$x]['longitude']."]";
    if ($x <= $num_of_rows - 2)
    {
        echo ",";
    }
  }

  echo "];";

}
vData();
?>

HTML script

 <!DOCTYPE html>

    <title>Tobago Street Map </title>
    <link rel="stylesheet" href="http://ajax.googleapis.com/ajax/libs/jqueryui/1.8/themes/black-tie/jquery-ui.css" type="text/css" />

    <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css" />


texts: GeoServer DescribeFeatureType response missing attribute for GML only

texts:
We have a GeoServer 2.18.2 instance serving data from PostGIS that seems to return a different structure for DescribeFeatureType requests depending on the outputFormat specified in the request.
One of our database tables which is published via GeoServer has 7 column names: [id, geom, code, name, enw, inset_no, category], nothing out of the ordinary. Making a WFS DescribeFeatureType to GeoServer for this and including outputFormat=application/json provides a complete description of the feature including all of the above attributes from the underlying table. However, if outputFormat is omitted so that GeoServer returns a GML3.2 response, then the name attribute is completely missing from the result.
There are no errors in the log from GeoServer when making these requests and any GetFeature requests against the same data source include the name attribute as expected. It is only DescribeFeatureType requests that seem to be affected, which causes a problem for some clients such as QGIS that use a DescribeFeatureType call to set up a layer. Prefixing the database column and reloading the feature type in GeoServer makes the attribute appear in responses but this feels like a workaround.
Is this expected behaviour from GeoServer or am I seeing a bug?


texts: How to tell which menu/toolbar a plugin gets added to?

texts:
I find a common issue with QGIS is installing the "Really Helpful Plugin You Really Needed", then spending 10-15 minutes in GUI trying to find it.

Is there somewhere, either some section in plugin documentation or in the log messages panel, that tells you the location within the GUI where the plugin was installed?


texts: Update EVLR header info for very large LAS files using laspy

texts:
I ran into another problem when applying my workaround to large LAS files.
The issue is noted at the github page for laspy here: https://github.com/grantbrown/laspy/issues/51
This builds on an earlier post at: Update EVLR header info for .LAS file using laspy

My method shown at the post above works for smaller LAS files (tested on files up to 800MB), but fails on larger files (fails on a file 6.2GB in size).

The issue seems to be the file byte offset integer value (where the EVLR content is to be placed) is too large for laspy (using version 1.2.5 with python version 2.7.12 64-bit) to convert to a C long value.

My error stack trace looks like this:

Traceback (most recent call last):

File "C:\Projects\LASPY_EVLR\writeEVLR.py", line 112, in <module>
  main(sys.argv[1:])
File "C:\Projects\LASPY_EVLR\writeEVLR.py", line 103, in main
outFile_v14.header.evlrs = old_evlrs
File "C:\Users\am\AppData\Local\Continuum\Anaconda2\envs\lidartools\lib\site-packages\laspy\header.py", line 1146, in set_evlrs
self.reader.set_evlrs(value)
File "C:\Users\am\AppData\Local\Continuum\Anaconda2\envs\lidartools\lib\site-packages\laspy\base.py", line 848, in set_evlrs
dat_part_1 = self.data_provider.fileref.read(old_offset)

OverflowError: Python int too large to convert to C long

Digging deeper, I found that the offset value for the beginning of the EVLR content is set to a value of 6,440,746,607 or about 6.44e9. I know this WILL fail to be converted to a double, but I'm not sure about C long types.

If this is a 32-bit unsigned long integer, then the conversion won't work... Is it possible to make this a 64-bit unsigned long integer? It looks like the max 32-bit unsigned long integer value is: 4,294,967,295.

This is a very large LAS file, though we have tons more of this size. This LAS file has 214691510 points.

I'm setting the "start_first_evlr" property to the end of the input file (when there are no EVLRs present) in the code from my previous post:

endofthefile = theDataProvider._mmap.size()
outFile_v14.header.start_first_evlr = endofthefile

Then when I set the updated EVLR content, the Overflow Error shows up:

outFile_v14.header.evlrs = old_evlrs

Do you have any workaround suggestions for me?


texts: Finding metadata for Fusion Table used in Google Earth Engine

texts:
I'm just starting to learn GEE, and I see a number of places use a Fusion Table that has polygons of world countries. It's used in the Google Earth Engine manual, tutorials, and questions on stackexchange. The Fusion Table is ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw and can (currently) be found here. But I can't find any reasonable metadata for it. I've tried:

getting the metadata within GEE:

// Country Fusion Table
var countries =  ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw');

// Get a list of all metadata properties.
var properties = countries.propertyNames();
print('Metadata properties: ', properties); // ee.List of metadata properties

// Get a specific metadata property.
var name = countries.get('name');
print(name); // ee.Number
var docid = countries.get('DocID');
print(docid);

But all that it returns is:

Metadata properties:
["name","DocID"]
   0: name
   1: DocID
Countries.csv
1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw

Going directly to the Fusion Table and going to file > about, which shows this:

Which still says nothing of where the info actually comes from or who created it.

So how do I find out more info about this particular dataset, and how do I evaluate the data quality of Fusion Tables more generally?


texts: deploying geoserver on linux or windows: performance implications?

texts:
I am deploying my first geoserver in production. I develop on a Macbook Pro retina, and my client uses Windows, so I deploy on Windows servers.

My first impression was that geoserver seemed a lot slower on the Windows server than my development machine, but I was able to improve this greatly using all the performance tips found all over the net. In short: install Native JAI, tune jvm, enable gwc. But still.

In my search for performance tips, I stumbled accross the following results: http://www.slideshare.net/gatewaygeomatics.com/wms-performance-shootout-2011

In this slideshow is shown that MapServer on linux is considerably faster than Mapserver on windows. I tried to look for a similar results or benchmark for Geoserver, but did not find it.

Aside of that I also wonder if the container has any impact on the performance of GeoServer (e.g. jetty/tomcat/jboss/...).


texts: Grouping or Clustering Points that are within 45 miles of each other using ArcGIS Desktop?

texts:
I have a project in ArcGIS (ArcMap 10.2) where I want to identify naturally occurring dense spots of points and then form groups in that dense area, where points are within 45 miles of each other. This would mean that not so dense areas will be left out and not become part of any group?

I don't have access to Point Density tool. I run Python code for Grouping Analysis but I get errors due to Python or GIS version mismatch.


texts: Viewing/editing attributed relationship data in Identify/Attributes dialog in ArcMap?

texts: Is there any extension or standard  option to view/edit attributed relationship records attributes using standard ArcMap attribute editor and/or Identify windows?


texts:
After opening attribute table, and right click the column header, I get a prompt to organise columns.

Expected: clicking and dragging a column on top of another would put it on top

Actual: column goes way under

Is there another way of organising columns in QGIS ?

This is happening on QGIS 3.4, 3.8 and even on 3.10

Edit: this is on mac version as people in the office use mac


texts: How to export a map to a single page pdf?

texts: I want to export my 900 x 800 mm map to a PDF file and all I'm getting is a file with 31 page of parts of my map. How do I create a single page - whole map pdf?


texts:
I'm writing a plugin using QGIS. I want to associate a shortcut to a procedure called "test".

The code i am using is the following :

self.shortcut = QShortcut(QKeySequence(Qt.Key_Q), self.iface.mainWindow())
self.shortcut.setContext(Qt.ApplicationShortcut)
self.shortcut.activated.connect(self.test)

For unknown reasons it seems to work erratically:

When I define a new key sequence (let's say "A") it is ok. I change the code, reload the plugin (with plugin reloader tool) the shortcut is working.
Then the shortcut seems to get unresponsive (I don't find the reason why..)

I choose another key sequence (let's say "Z"), it is working.

When I try to switch it back to a previous setting (back to "A" for example) it is not working either ...

It seems that only new key sequence choices are working (for a while only ...)

What could I be doing wrong?

Is it an issue that could occur only during design mode of the plug-in?


texts:
https://www.dropbox.com/s/disaisc1wrv1mbi/Form470_2015-02-11-13-10-13.xls?dl=0

I downloaded that from the USAC (people that manage USF, E-Rate, etc.) web site. I've tried many combinations of importing to Access, linking ArcCatalog to Access, ensuring headers don't start with a number, changing spaces to underscores in headers, etc. I always end up with:

Start Time: Wed Feb 11 12:32:14 2015
ERROR 999999: Error executing function.
Failed to execute (TableToTable).
Failed at Wed Feb 11 12:32:14 2015 (Elapsed Time: 0.00 seconds)

I read a few posts on different sites about the issue, but it was always with people using their own Python scripts. I'm just using the standard ArcCatalog import function.

Ideas?


texts: Blending overlapping raster mapsheets with QGIS and GIMP

texts:
How can overlapping raster map sheets be overlaid in QGIS, and manually blended together in GIMP?

I have an overlapping set of hand-drawn raster maps that contain varying levels of detail and would would like to combine these into a single large map using QGIS and GIMP. Since the information in the overlap regions is sometimes complementary, I would like to be able to choose manually which pixels to take from which map sheet.

I was successful at georeferencing them onto OpenStreetMaps, and would now like to blend the overlaps together (which I cannot seem to work out how to do in QGIS). In GIMP, using one layer for each map of a given two map sheets, I can delete the bits in the top layer which I want to use from the lower one, make these regions transparent using a layer mask, and then merge down.

How do I transfer my aligned map sheets as individual layers into GIMP for editing and then back to QGIS?

I am new to QGIS and have played around with this for a while now but can't quite figure it out. Maybe I'm overlooking something or I am on the wrong track altogether.


texts: Arcgis apply lyr file to multiple layers in workspace

texts:
I've adapted some code from Apply same symbology to many raster layers in ArcMap? and the ESRI noticeboards to batch apply a .lyr file to several hundred .TIF files:

import arcpy
#Specifies mxd
mxd = arcpy.mapping.MapDocument(r"PATHTOMXD\ModelVisuals.mxd")

#Specifies layer file with the symbology you want to apply to the other raster layers
sourceLayer = arcpy.mapping.Layer(r"PATHTOLAYERFILE\FloodDepths.lyr")

#Loops through layers in the mxd
for lyr in arcpy.mapping.ListLayers(mxd,""):
    arcpy.ApplySymbologyFromLayer_management(lyr, sourceLayer)

The code runs fine with no errors and it loops through the rasters, locking them one by one, but it doesn't apply the symbology layer or turn them on.

I am very rusty with python for ArcGIS.  What am I doing wrong or is there a better way of doing it?

####

Update:

This code (from this question) works perfectly:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.activeDataFrame
groupName = "Group"
lyrFile = "XXXXXXX"
allLayers = arcpy.mapping.ListLayers(mxd, groupName, df)[0]
for layer in allLayers:
     arcpy.ApplySymbologyFromLayer_management(layer, lyrFile)
     layer.visible = "True"
arcpy.RefreshActiveView()

What is the purpose of the "[0]" at the end of the allLayers line?  Although it is in the code snippets here it isn't explained.


texts: Classifying neighbors AND neighbors of neighbors using ArcGIS Desktop?

texts:
In ArcGIS 10.2, I'm trying to identify the neighbors AND neighbors of neighbors for counties in the United States.

Although I can identify which counties are neighbors to each county using either a distance table or the Polygon neighbors tool, I'm having trouble identifying the neighbors of neighbors while also excluding those neighbors of neighbors that are neighbors of the county.

Here is a picture of what I'm trying to identify. I'd like to identify the light blue features for each red feature. The dark blue counties are the neighbors of the red county; I want to identify the light blue one (i.e., neighbors OF neighbors of a feature)...


texts: Export Attachments using Python script to modify file name [ArcGIS Desktop 10.8.1]

texts:
I'm trying to export photo attachments from a railroad crossings dataset taken with ArcGIS Collector, and would like to customize the naming convention so that the file name contains the crossing number. I stumbled across a script from @Midavalo that accomplishes this, however the files exported are no longer JPEG format. I have almost no knowledge of Python.
Feature table:

Attachment table:

Modified script:
import arcpy, os
from collections import defaultdict

inFC = r'C:\Users\Desktop\scratch\CustomExportAttachments\CopyFeatures.gdb\Crossings'
inTable = r'C:\Users\Desktop\scratch\CustomExportAttachments\CopyFeatures.gdb\Crossings__ATTACH'
fileLocation = r'C:\Users\Desktop\scratch\CustomExportAttachments\Attachments'

# Get dictionary of ObjectID and associated field value
myFeatures = dict()
with arcpy.da.SearchCursor(inFC, ['OID@', 'CROSSING']) as cursor:
    for row in cursor:
        myFeatures[row[0]] = row[1]

# Create dictionary to count usage of the field value (to increment files)
valueUsage = defaultdict(int)

# Loop through attachments, incrementing field value usage, and using that
# increment value in the filename
with arcpy.da.SearchCursor(inTable, ['ATTACHMENTID', 'ATT_NAME', 'REL_OBJECTID']) as cursor:
    for row in cursor:
        if row[3] in myFeatures:
            attachment = row[0]
            fieldValue = myFeatures[row[3]] # Value of specified field in feature class
            valueUsage[fieldValue] += 1 # Increment value
            filename = "{0}_{1}".format(fieldValue, valueUsage[fieldValue]) # filename = FieldValue_1
            output = os.path.join(fileLocation, filename) # Create output filepath
            open(output, 'wb').write(attachment.tobytes()) # Output attachment to file

Output folder:

Why are the photo attachments exporting as "File" instead of JPEG, and how should I change the script in order to make this happen?


texts: Split a polygon using SqlServer geometry type

texts:
I have a polygon geometry in SqlServer and i want to split that single polygon into two features via a line that i draw through it.
How can i do that with sql?
I don't want to lose ANY portion of the original polygon. So buffering the line and doing a difference is not going to get it done. Unless there is a way i can "put back" the lost edges from the diff.

SqlGeometry splitLineGeom = SqlGeometry.STGeomFromText(new SqlChars(splitLineWKT.ToCharArray()), Srid);
// buffer line and get the difference from the input
SqlGeometry bufferedSplitLine = splitLineGeom.STBuffer(0.001);
SqlGeometry resultGeom = inputGeom.STDifference(bufferedSplitLine);


texts: Directory.GetFiles vs. ListDatasets - Which is the tool for the job?

texts:
I'm currently writing an arcobjects app that performs an arbitrary set of geoprocessing tasks on a group of shapefiles and/or dbf tables. My current implementation uses Directory.GetFiles to build an array of file paths which I feed to various GP execution functions. I'm wondering if ListDatasets can perform the same task. I need to search an entire directory tree starting at a given root folder for all files which match a given search string. I am not certain if ListDatasets will do this (I have a devil of a time deciphering the SDK documentation).

In short, which would prove more efficient and effective and recursively searching a directory tree for .shp and .dbf files, Directory.GetFiles, GP.ListDatasets, or some heretofore unmentioned method?


texts: Error: Input [, ] is the wrong shape for a LineString

texts:
I'm trying to create a linestring out of a dataframe of latitudes and longitudes.

geometry = [Point(xy) for xy in zip(linkdf.lon, linkdf.lat)]
linkgdf = GeoDataFrame(linkdf, geometry=geometry)
linkdf_line = linkgdf.groupby(['linkId'])['geometry'].apply(lambda x: x.tolist())
linkdf_line = linkdf_line.apply(lambda x: LineString(x))
linkgdf_line = GeoDataFrame(linkdf_line, geometry='geometry')

Where linkdf is just a dataframe of latitudes and longitudes:

Unfortunately, I keep getting this error:

ValueError: Input [<shapely.geometry.point.Point object at 0x000002AB98131AC8>, <shapely.geometry.point.Point object at 0x000002AB98131198>] is the wrong shape for a LineString

What is really confusing me is that this seems to work in a jupyter notebook, but not running my python script in terminal.


texts: How to find which direction to flow between two points in QGIS Python?

texts:
I want to find which way direction to go based on two points and line segment.
I have two layer line layer and point layer.user selects the two points in lines layer.
I have to find which direction to go from first point to second point based on line segments to collect intermediate points.

All the lines in line layer is breaked every point to point.I want only one direction .

I have already tried create the buffer for first point to get inital two lines and loop through direction left or right until second point found. It takes more time.

How can I do that?

.


texts: Visual Studio 2019 intellisense for arcpy not working in code window

texts: Using Visual Studio 2019 intellisense works for the immediate window but not the coding window. How can I activate intellisense for coding window?


texts: Exclusive Economic Zone (EEZ) and Regular Boundaries one shapefile

texts: Does anyone have an idea if there is a shapefile of the merged EEZ boundaries and regular administrative boundaries? I need this data for entire Africa. I have downloaded the EEZ boundaries from http://www.marineregions.org/downloads.php but they don't match with any of the administrative datasets that I have located that show the African continent boundaries. There are gaps and overlaps. I can do merge or similar but since it is for an entire continent, it would take a lot of time? I am surprised not being able to find the merged boundaries online. Any ideas?


texts: Correlation/association between two point layers

texts:
I’d like to know if there is a way of obtaining a correlation/association test between two point layers in QGIS (or another open source software if not possible in QGIS)?

I have two point layers- one layer is the location of individuals (long and lat coordinates) who have contracted an infection. The second layer is the long and lat coordinates of important buildings (e.g. schools). I’d like to know if the location of infected individuals is correlated/associated with the location of important buildings.


texts: Filter the Feature Layer parameter in Python Script Tool

texts:
I have a Python script tool where I have a parameter of type Feature Layer because I want to have the already opened layers as a dropdown list.

I was wondering if there is a way to filter this Feature Layer parameter dropdown list in order to show Polygon Shapefiles only?


texts: Is it possible to convert QGIS 3 measure tool output to polyline?

texts: I am wondering if it is possible to convert the results of the QGIS (3.16.0-Hannover) measure tool to a polyline (or shapefile). I need to make a series of linear measurements and would like to keep a record of the geometry of what I measured. When I try to paste the results of the "Copy All" from the measure tool I get an empty attribute table however the number of blank features shown is the same as the number of measurement segments so, it feels like the info is there... somehow....


texts:
I have created an azimuth orthographic projection map (basically a globe, I've used this custom CRS: +proj=ortho +lat_0=51.470129 +lon_0=-0.452751 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs with QGIS 3.4.4. The labels are flat, as you would expect for readability purpose. Does anybody know if it's possible to have the labels following the earth surface curvature instead, like if it was an actual picture of a physical globe?

If it's not possible with QGIS, how about other software programs?


texts: Extracting mean elevation values from DEM based on multi-polygon shapefile?

texts:
I'm trying to get the overall average elevation for different areas delineated within a polygon.
Originally, I converted the Digital Elevation Raster (DEM) to a polygon and tried to merge the attribute tables. However, the attribute table for the DEM doesn't actually contain the elevation data.
So,then I converted the DEM into points and extracted the elevation values for each of these points using the 'Extract multivalues to points'....which saved the elevation to the point file.
Now having done that, I don't know how I can merge the point file with my polgon file to get the mean elevation for each of the individual polygons. If I mask the point file using the polygon, I'll just get loads of points per polygon and I want a single value.

The picture shows the individual areas I want the elevation for.
Can anyone advise?


texts: Geoserver: Call custom OWS function from inside SLD

texts:
I've a PostGIS table published  as a layer, I need to style this layer using SLD.
I need to call a custom OWS function (which I wrote) from inside the SLD. It should look something like:

<?xml version="1.0" encoding="ISO-8859-1"?>
<StyledLayerDescriptor version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/sld http://schemas.opengis.net/sld/1.0.0/StyledLayerDescriptor.xsd"
  xmlns="http://www.opengis.net/sld" xmlns:ogc="http://www.opengis.net/ogc"
  xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">

  <NamedLayer>
    <Name>test_fn_call</Name>
    <UserStyle>
      <Title>dark yellow square point style</Title>
      <FeatureTypeStyle>
        <Rule>
          <Title>dark yellow point</Title>
          <PointSymbolizer>
            <Graphic>
              <Mark>
                <WellKnownName>square</WellKnownName>
                <Fill>
                  <CssParameter name="fill">#99cc00</CssParameter>
                </Fill>
              </Mark>
              <Size>
              <ogc:Function name="gs:calcSize">
                <ogc:Function name="parameter">
                    <!-- **** What should I write here? **** -->
                </ogc:Function>
              </ogc:Function>
              </Size>
            </Graphic>
          </PointSymbolizer>
        </Rule>

      </FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>

calcSize is a function to calculate point size depending on some calculations on the feature's fields.

My question is what should I write as a parameter of this function, actually I need to path the current feature as a parameter, so what should I write in the SLD as a parameter?

Also in the java side, what should the parameter type be?


texts: Failure to connect QGIS with GeoServer WMS?

texts:
I'm trying to add a WMS to QGIS from GeoServer from localhost but after I paste in the URL and go to connect I get the following error:

Failed to download capabilities:
Download of capabilities failed: Error downloading http://localhost:8080/geoserver/geog585?SERVICE=WMS&REQUEST=GetCapabilities - server replied: Not Found

I'm pretty sure I typed in the right parameters but I'm a total webmap nube so I could be off:


texts:
I am using the Google Earth plugin for QGIS.  I have the following code that returns the DEM for Oman.
    import ee;
    image = ee.Image("USGS/SRTMGL1_003");
    Oman = ee.FeatureCollection("FAO/GAUL_SIMPLIFIED_500m/2015/level0");
    filter = ee.Filter.eq('ADM0_NAME','Oman');
    OmanFiltered = Oman.filter(filter);
    Map.addLayer(image.clip(OmanFiltered),{'palette': ['blue', 'red'], 'min': 0, 'max': 5000},'OmanDem',True);

The problem arises in that this resultant map has no extents being based on what is called a 'global composite' in GEE.  I cannot therefore Export->Save As the layer as a GeoTIFF with an error saying 'Attempt to create 0x0 dataset is illegal'.  I guess that is because it has no extents.
I tried to draw a vector polygon in QGIS as a surround of the GEE supplied image above and then do a Raster->Extraction->Clip Raster by Mask Layer using the surrounding polygon but then got an error
 ERROR 4: `type=xyz&url=https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/maps/c65a4d164f5547eb7fb9f7fe1c1b24df-c4621e8e2576f3167a45a4a52e46ec04/tiles/{z}/{x}/{y}' does not exist in the file system, and is not recognized as a supported dataset name.

Does anyone have a solution to extracting the raster files from GEE in QGIS that allows them to be saved as a GeoTiff for further processing in QGIS.


texts: Adding QGIS feature, that shows coordinates and scale, into a custom application

texts:
I'm writing a standalone application using PyQGIS and I'd like to add the QGIS feature that shows coordinates and scale as the cursor hovers over the map.

I need to know this mainly in case of raster maps.

Can someone guide me?


texts: QGIS: Burning stream network (catchment area)

texts:
I need help from stackexchange because of a DEM correction problem...

I have a DEM with a 5 meters resolution where one of the main river is cut by a large highway. In order to delete that topology error I'm looking for a method for burning the river in the DEM.

After looking at tutorials I've tried "r.carve" (additionnal stream depth = 10 meters) but even thought it happened to have corrected my DEM, the flow accumulation in this area is very discontinuous...

What is the problem here?


texts: gdal_translate syntax windows command line

texts:
I am trying to extract a band from a MODIS hdf file on windows with gdal_translate.

Here's my prompt:

gdal_translate -of GTiff "HDF4_EOS:EOS_GRID:\"MOD13A2.A2009001.h17v08.005.2009020092120.hdf\":MODIS_Grid_16DAY_1km_VI:1km 16 days NDVI" NDVI.tif

However I get:

ERROR 4: `'HDF4_EOS:EOS_GRID:MOD11A2.A2009001.h17v07.005.2009014044746.hdf:MODIS
_Grid_8Day_1km_LST:LST_Day_1km'' does not exist in the file system,
and is not recognised as a supported dataset name.

GDALOpen failed - 4
`'HDF4_EOS:EOS_GRID:MOD11A2.A2009001.h17v07.005.2009014044746.hdf:MODIS_Grid_8Da
y_1km_LST:LST_Day_1km'' does not exist in the file system,
and is not recognised as a supported dataset name.

I imagine it's something to do with my syntax being a bit 'iffy' though it seems perfectly fine to me

Help much appreciated


texts: How to write point's attribute to polygon table using pyqgis?

texts:
test2 is point data and test is polygon data
I will write test2's attribute to test table i want to make it

 i want

 date|Num|add|결과|x|y|BUFF_DIST|contain
 20170103|640272|강원도 인제군...|정좌표|318102|501346|30000.0000|86,91,235,2485...|

I want to use pyqgis and I wrote code but I don't know how to write point's attirbute to polygon table

my code is

#-*- coding: utf-8 -*-
from qgis.utils import iface
from PyQt4.QtCore import QVariant
NF = "FRMHS_NO"
f = iface.addVectorLayer("D:/2018/nonmoon/test/test.shp","test","ogr")
g = iface.addVectorLayer("D:/2018/nonmoon/test/test2.shp","test2","ogr")
f.startEditing()
res = f.dataProvider().addAttributes([QgsField("contain",QVariant.String)])
f.updateFields()

I have a image. first img is map
second img is test(polygon data)'s table
last img is test2(point data)'s table


texts: Using custom map tiles with Google Map API V2 for Android.?

texts:
I'm looking for a way to use custom map tiles with the Google Map API V2 for Android.

I'm writing an application that will create it's own maps in real time with data that are comming from a robot.

The application need to show to the operator this map. The operator need to interact with this map, letting way points, etc.

I want to use the GoogleMap engine in order to do the same that this page :

http://cdn.mikecouturier.com/blog.mikecouturier.com/tilesgenerator/index.html

The problem is that he used the Javascript API when I want to use the Android API

Is there any way to use Custom tiles map on android with the Google Maps Engine?

I'm already looking at how to use ArcGIS but, i prefer using an API without paying for a license.


texts: Semi-Automatic Classification Plugin Image preprocessing Error

texts:
I'm using QGIS 3.10 for MacOS with Semi-Automatic Classification Plugin 6.4.7
Currently, I'm trying to perform initial image pre-processing of one of the Sentinel-2 images.

However, when I click Run I got the following error

Has anybody encountered such or similar problem?
Unfortunately, Google didn't help me this time.


texts: Using "If Value Is" tool in ModelBuilder of ArcGIS Pro?

texts:
I want to compare an input number to a constant number (<=), but I can't find the tool to do it in ModelBuilder. I have seen in an ArcGIS help website that a function named "if value is" exists but I don't know how to find it.

Do anyone know how to do it?

https://pro.arcgis.com/en/pro-app/tool-reference/modelbuilder-toolbox/examples-of-using-logical-tools-in-modelbuilder.htm


texts: How to improve precision of geographic coordinates using ogr2ogr to translate from shapefile to KML?

texts:
I made a shapefile in QGIS (ver 3.16) project with EPSG:4326. The shapefile is for a simple polygon and I want to display the polygon on Google Earth so I need to convert it to KML format. The default CRS in Google Earth is also EPSG:4326.
I am using ogr2ogr to perform the conversion:
ogr2ogr -t_srs EPSG:4326 -f KML test_polygon.kml test_polygon.shp

But visually the position of the KML in Google Earth is offset by about 0.10 second to the west (about 10 feet), and about 0.01 second to the south from the position shown in QGIS. See images below (ogr2ogr does not retain styling, so the KML only shows the outline of the shape, but you can still see the position offset).
Shapefile position shown in QGIS:

KML position shown in Google Earth:

Is it possible to improve the precision of the conversion to KML such that the position shown on Google Earth is the same as shown in QGIS? Or perhaps I am going about this the wrong way?


texts: Pass Queried Feature Set to geoprocessor as a parameter

texts:
I am new to geoprocessing services and tried to publish my first.  I created a simple model with one parameter, a featureset, and I am intersecting it with another layer which is specified in the model and uploaded to the server.

If I pass the required parameter,as below, it never executes.  If I pass an empty parameter object, var params = {}, it runs.

    var app = {};
    app.qtCounties = new QueryTask("http://myserver/arcgis/rest/services/PA_Counties/MapServer/0");
    app.qCounties = new Query();
    app.qCounties.returnGeometry = true;
    app.qCounties.outFields = ["*"];

    paCounties.on("click", getData);function getData(evt) {
            console.log("In GetData The Event graphic is " + evt.Graphic + "Point is ", evt.mapPoint.x, evt.mapPoint.y)
            var county, promises, qGeom, point, pxWidth, padding;

            point = evt.mapPoint;
            pxWidth = map.extent.getWidth() / map.width;
            padding = 2 * pxWidth;
            qGeom = new Extent({
                "xmin": point.x - padding,
                "ymin": point.y - padding,
                "xmax": point.x + padding,
                "ymax": point.y + padding,
                "spatialReference": point.spatialReference
            });

            app.qCounties.geometry = qGeom;

            county = app.qtCounties.execute(app.qCounties);
            console.log("deferreds: ", county);
            promises = new all([county]);
            promises.then(handleQueryResults);
            console.log("created d list");

            function handleQueryResults(results) {
                console.log("queries finished: ", results);

                gp = new Geoprocessor("http://myserver/arcgis/rest/services/TestIntersectModel/GPServer/TestIntersectModel");
                gp.setOutSpatialReference({ wkid: 4326 });

                var params = { "Input_Polygon": results[0].features };
                gp.execute(params, completed);
                function completed(res) {
                    console.log("Completed GeoProcssing", res)
                }

                console.log("Params ", params)}

When the query runs I get a valid object in results which is set as the lone parameter. Apparently I am not passing the parameter as a featureset but not sure how else to pass it?  Any help would be appreciated.

I published my service using this example.  Though, I did not publish it as a map service, the last step, and I am trying to access it via javascript api instead of the server web page method.  The model works fine when run in desktop and providing the parameter manually.


texts:
I have a points layer with bus stop locations and a line layer which represents the road network. I want to find the shortest path for each bus route based on the sequence of bus stops and save the path as a new layer or extract the paths into a csv file.
I have ensured that the bus stops fall exactly on the road network.
How do I do this in QGIS?


texts:
I'm starting to explore ArcGIS Online as my first introduction to web mapping (perhaps not the best choice of cloud services but thought I'd give it a shot). I was intrigued by their discussion of "story maps" and wanted to create one that overlaid raster imagery that I had on my system. However, I can't figure out how to add raster data. It seems to only accept .shp, .csv etc. I can add rasters if they're converted to layers, but can't seem to actually overlay them on the map. In looking through the help it seems the option is using a map service, the more I dig into this the more complicated it seems.

Is it not possible to simply overlay raster data as one would do in ArcMap Desktop?


texts: R spatial: st_overlaps to subset also fully covered polygons?

texts:
I want to subset all polygons in fc1 completely or partially overlapping with fc2. But, the st_overlaps subsets only polygons that overlaps with the border of fc2, but skips the fully covered polygons. I.e. here: polygons in center (grey) are not subsetted, only red polygons are returned as overlapped:

I wonder if there is a way how to subset all - partially and fully covered - polygons at once? Is really the only option to subset polygons twice: by st_covers and st_overlaps and then combine them? Would not this approach lead to duplicated polygons?

My dummy example:

# Load data
shp = system.file("shape/nc.shp", package="sf")
my.sf <- st_read(shp, quiet = TRUE)

# Convert crs to projected system to make buffer
my.sf.web<- st_transform(my.sf, 3857)

# Subset one polygon
i = 10
one  = my.sf.web[i, ]

# Create buffer
buff50 = st_buffer(one, 50000)

# Check which polygons overlaps with my buffer
out.overlap = st_overlaps(buff50, my.sf.web)

# Subset the polygons that overlaps with the buffer
nbrs.buff.over  <- my.sf.web[st_overlaps(buff50,my.sf.web)[[1]],]
nbrs.buff.cover <- my.sf.web[st_covers(buff50,my.sf.web)[[1]],]

# Merge both spatial subsets:  # !!! Is there easier and less error prone solution??
out <- cbind(nbrs.buff.over, nbrs.buff.cover)

I expect that my output will contain all - completely and fully covered polygons.


texts: Calculating Edge Percentage using ArcMap?

texts:
I have some NCLD data of an area I've clipped out and there is a count of each habitat class (open water, deciduous forest, etc.).

I'm trying to find a way to calculate the edge amount of these habitat classes in ArcGIS 10.6.

Here's a screenshot of what I'm talking about:


texts: OGR GeometryReference Intersects() method with OSGB projection causing crash

texts:
I'm trying to use the GeometryReference.Intersects() method to see if a feature within a grid intersects with another polygon. I have tried it in the QGIS python console (windows) and directly on linux, and it crashes in both cases.

EDIT: I have tried the same code with the data reprojected to UTM, which is working. So possibly it is a problem with using the Intersect method the the OSGB projection (EPSG: 27700)

If possible I would like to know what is going wrong, or alternatively if there is another method of checking if two features intersect within python.

Here is the code I am trying:

from osgeo import ogr

driver = ogr.GetDriverByName("ESRI Shapefile")

grid_ds = driver.Open("test_grid.shp")

poly_ds = driver.Open("test_poly.shp")

grid_layer = grid_ds.GetLayer()
poly_layer = poly_ds.GetLayer()

poly_geom = poly_layer[0].GetGeometryRef()
grid_geom0 = grid_layer[0].GetGeometryref()

grid_geom0.Intersects(poly_geom)

Here is a link to some example datasets


texts: Obtaining length of river upstream of points up to next point using ArcGIS Desktop?

texts:
I am carrying out a study which aim is to come up with a list of priority obstacles to fish migration to be removed based on those where more extra habitat would be created after the modification.

I have a river network and the location of some dams, and I would like to know the length of the river available upstream of each dam up to the next obstacle. I need to obtain the total length not only of the main river, but also the tributaries for each dam.

The dams are the blue points and I would like to obtain the total length of the river reaches within the boxes. Some of the points aren't on the river itself (precision of the dam location might not be great) so a buffer of 100m should be applied to link them to the nearest river. I am using the OS Open Rivers Network, but I also have the OS Water Network if that helps at all.

I have tried the Network analyst but it hasn't quite worked. The example below shows some of the issues:

I would manually have to identify start and end points for each of
the obstacles, and there are over 5000 of them, so I would like to
automate the process.
It measures the distance between nodes rather than the points I have specified, so it includes some river
downstream from the first dam, and it doesn't go all the way to the
upstream obstacles.
Some reaches are completely missing.

I have tried to add the obstructions into the nodes info as "Pseudo" in category but when I created the geometric network, but it returned 1000s of errors. Some of them are type 11 (feature's geometry has multiple parts), but most of them are type 16 (junction isn't connected to any other edge feature). I don't get any of these errors when I use the original Hydronode.

Is there another way to include the obstructions into the Hydronodes that doesn't create?

I have had a look at the geometric network editor, but I would have to change every point individually or could it be done in bulk?

Could it be possible to create a network using the obstructions only as "Pseudo" or does it need junctions, outlets and sources?

I also have a polygon river layer and I would like to know not only the river length upstream of a dam, but also the river area, which can be translated into habitat available to fish. This is because it would be more beneficial to remove a dam in a river 10m wide than only 2m wide. However, this last point isn't as critical as obtaining the river length at least.


texts: Example of how to display a Cartodb layer (using the layergroupid) after initializing a named map?

texts:
I have successfully created a named map. I can initialize the map (as seen here), and I do get the correct response. Now, to my problem.

How do I show this layer (using the layergroupid I got as a response) on a map with CartoDB.js?

When I look at the documentation I see:

  You can then use the layergroupid for fetching tiles and grids as you would normally (see anonymous map section)

In the anonymous section, I see the following

  The tiles can be accessed using:
  http://{account}.cartodb.com/api/v1/map/:layergroupid/{z}/{x}/{y}.png

How is this done on the client side? Are there any examples out there where you start with a blank map, initialize a named map and add that layergroupid to the map? I have not managed to find such a thing.


texts:
I got a dwg-file which I imported in QGIS, it is displayed allright.
There was a prj-file which stated that that the CRS should be EPSG 31258, which I used for importing the file. However it does no match at all, it it isn't even in the area of OSM.

I checked the file in TrueView, there it shows "WCS".
There coordinates displayed in Qgis and TrueView are about: 89708467,340439859.
All common CRSs (at least the one I  know) for the areas were checked.

Is there a way to determine or set the correct CRS?


texts: How to show label of a new feature layer in OpenLayers?

texts:
I have created a dynamic new feature layer along with a marker using OpenLayers with popups, but, I want to show labels for each new feature created...

Below is the code I used...

//To add new Vector Layer
var style = new OpenLayers.Style({
    fillColor: "#ffcc66",
    strokeColor: "#ff9933",
    strokeWidth: 2,
    label: "${type,area}",
    fontColor: "#333333",
    fontFamily: "sans-serif",
    fontWeight: "bold"
});
var type='MAHAMAYATALA';
var vectorLayer = new OpenLayers.Layer.Vector("SupplyOff",{styleMap: new OpenLayers.StyleMap(style)});
//  For Consumer Query
$qry_xy="select * from ci_master@ltdvlp  where  x is not null  AND GIS_SOURCE_CODE LIKE '01%' AND ROWNUM < 5";
$parsed =oci_parse($db_conn_prod,$qry_xy);
oci_execute($parsed);
$nrows3=oci_fetch_all($parsed,$results3);
echo $nrows3;
for($i=0;$i<$nrows3;$i++){
$x=$results3["X"][$i];
$y=$results3["Y"][$i];
$consno=$results3["CONS_NUM"][$i];
$source=$results3["GIS_SOURCE_CODE"][$i]; ?>
        var x = '<?php echo $x;  ?>';
        var y ='<?php echo $y;  ?>';
        var cons ='<?php echo $consno; ?>';
        var type = '<?php echo $source; ?>';
         // Define markers as "features" of the vector layer:
var feature = new OpenLayers.Feature.Vector(
        new OpenLayers.Geometry.Point( x, y ),
        {description:'cons: '+cons} ,
        {externalGraphic: 'img/marker-gold.png', graphicHeight: 25, graphicWidth: 21, graphicXOffset:-12, graphicYOffset:-25  }
    );

vectorLayer.addFeatures(feature);
map.addLayer(vectorLayer);

// for Pop ups Window
var controls = {
      selector: new OpenLayers.Control.SelectFeature(vectorLayer, { onSelect: createPopup, onUnselect: destroyPopup })
    };

function createPopup(feature) {
  feature.popup = new OpenLayers.Popup.FramedCloud("pop",
      feature.geometry.getBounds().getCenterLonLat(),
      null,
      '
'+feature.attributes.description+'
',
      null,
      true,
      function() { controls['selector'].unselectAll(); }
  );
  //feature.popup.closeOnMove = true;
  map.addPopup(feature.popup);
}

function destroyPopup(feature) {
  feature.popup.destroy();
  feature.popup = null;
}

map.addControl(controls['selector']);
controls['selector'].activate();


texts:
I was wondering if someone could help me out on how to calculate runoff volume within an urban area ?

I have already produced different outputs:

Filled DEM
Catchment area
Channel network
flow direction
junctions
watershed basins
drainage basin
maximum flow path lentgh
channel network
upslope area

And I have also ranked the different surfaces according to an imperviousness gradient.

Now I would like to identify the area that would be the most vulnerable to a flooding. Unfortunately I have been unsuccessful so far.

I took a look over r.sim water that seems to be a promising tools for that but I'm facing issues with it !


texts: Problem aligning points near North Pole using Web Mercator with QGIS and Google Street View

texts:
This is my first time using QGIS 2.18.

I created a project with EPSG:3857.

Added OpenLayers Plug-in to use Google Street View as a map for the points and data I want to work on.

I love self-learning and learn-by-doing, I watched some vids and got an idea on how this program work.

I thought I'd make a sample to test if what I understood was right or not. I chose a point on Google Maps, right-clicked and chose "what's there" to get the coordinates. it is written as xxxxx.xxxxN, yyyyy.yyyyyE.

I added this in an excel file with 4 columns: ID, Name, Lat, Long. Saved as CSV. Opened it with notepad, saved as txt.

I then added a delimited text layer with the same ESPG:3857. The point in the file is almost at North Pole.

I tried to use some web-converters to know where the problem is but the best thing I got was like 400 km from the original location of my point.

My objective is to find a way to get coordinates from Google Maps, store in a text delimited file, add a text delimited layer and get the points where they should.

Can you help?


texts: Arithmetic on rasters with variable names in Map Algebra/Spatial Analyst?

texts:
I am a Python novice trying to figure out how to perform calculations on temporary rasters with variable names.

My attempt at a minimum working example is below, but I would eventually like to do this with thousands of larger grids, so I expect that I will need to keep everything temporary until the very end to increase efficiency.

The code below fails on the last line. The error message is:

  RuntimeError: ERROR 000732: Input Raster: Dataset bin2 does not exist
  or is not supported.

If I replace the last line with the permanent filenames:

tmpout = Minus("D:/tmp/tmp" + `year`,"D:/tmp/tmp" + `yearold`)

then it runs fine. I am running IDLE in ArcGIS 10.0 with Python 2.6.5 (r265:79096, Mar 19 2010, 21:48:26) [MSC v.1500 32 bit (Intel)] on win32.

import arcgisscripting
import arcpy
from arcpy import env
from arcpy.sa import *

arcpy.env.extent = arcpy.Extent(0,0,2,2)

arcpy.CreateRandomRaster_management("D:/tmp/","tmp1","","",1)
arcpy.CreateRandomRaster_management("D:/tmp/","tmp2","","",1)
bin1 = arcpy.Raster("D:/tmp/tmp1")
bin2 = arcpy.Raster("D:/tmp/tmp2")

for year in range(2, 3, 1) :
    yearold = year - 1
    binyear = "bin" + `year`
    binyearold = "bin" + `yearold`
    print "fails in next line"
    tmpout = Minus(binyear,binyearold)

I think the problem is that bin1 and bin2 are being interpreted as strings as opposed to pointer to temporary rasters, but I can't figure out how to make my intent understood. Replacing binyear with Raster(binyear), for example, did not help.


texts: How to save raster with projection into Python

texts:
I applied a raster processing that is already in a known coordinate system with this code,

import cv2
import numpy as np
from matplotlib import pyplot as src

img = cv2.imread(r'E:\2_PROJETS_DISK_E\threshold\621.tif',0)

# Otsu's thresholding after Gaussian filtering
blur = cv2.GaussianBlur(img,(5,5),0)
ret1,th1 = cv2.threshold(blur,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)

for i in xrange(1):
 plt.imshow(images[i*3+2],'gray')
 plt.xticks([]), plt.yticks([])
plt.show()

cv2.imwrite( r'E:\2_PROJETS_DISK_E\threshold\621-1.tif',th1);

but when I save my raster at the end of this script, I get a non-georeferenced TIFF raster

How can I keep the same coordinate system of the initial raster (without tranforming the output raster into local coordinates)

Since I am new to Python, and I have no knowledge of Python, i would like someone to correct me my script please.


texts: Importing, moving, rotating and scaling jpg images onto OSM maps in QGIS

texts:
I'm a master planner and I want to use QGIS for, among others, importing images (jpg, png, etc without any georeferences). These images will mostly be screenshots of existing structures in Google Earth.

Now when I import such an image the scale and location are far off to what and where I want it. Does anyone know a good method to insert the image on the wanted location on the OSM base layer, translate, rotate and proportionally scale it?


texts: InsertCursor on a event table (linked to linear referencing) not displayed directly on the map

texts:
I have a Python script (Python add-in) that creates new entries (with InsertCursor) in an event table that was used for linear referencing (through Display Route Events). The table is stored in a file GDB.

My first problem is linked to this question:

Refreshing open attribute table in ArcMap using ArcPy rather than ArcObjects?

the table is not reflecting the edition, except if I close the table tab and relaunch it, or if I sort an attribute of the table. I have tried all solutions exposed in this question. But that’s OK, I can deal with a manual refresh.

My main problem is that the entries I have created are “dead datas”: they are not represented in the map. My only solution is either to close the mxd and relaunch it, or edit them again, and then the segment created is shown.

I have tried to edit the layer (not the table), to update the entries directly after their creation (in the same script), to recreate a new event layer (display route events) and to change to "da" cursors, but the problem remained.


texts: PyQGIS Turn on panel not working once inside class (QDockWidget)

texts:
I am relatively new to PyQGIS and struggling to convert my simple code into a plugin class.
The code works when run from the Python window, but when I add it to def run(self): I get an error, which I am sure is something simple, but I just can't get it.
I am trying to get the plugin to display the GPSInformation window.
def run(self):
    for x in iface.mainWindow().findChildren(QDockWidget):
        if x.objectName() == 'GPSInformation':
            x.setVisible(False)

The error says NameError: name 'QDockWidget' is not defined
Thanks.


texts: Exporting only certain fields from a PostGIS table using QGIS Python

texts:
I'm trying to export tables from a PostGIS database into MapInfo TAB files, but only certain rows and columns. I have my rows exporting right, using a query, but I'm stuck on just getting the columns I'm interested in. I thought the

request.setSubsetOfAttributes(["obj","name","railway"],vlayer.pendingFields())

below would do the trick, but I'm still getting all of the original fields in my output.

    if table == "line_table":
        fields = vlayer.pendingFields()

        railway = ["Railway", r""""railway" != ''"""]
        request = QgsFeatureRequest()
        request.setFilterExpression(railway[1])
                   request.setSubsetOfAttributes(["obj","name","railway"],vlayer.pendingFields())

        selectedFeatures = vlayer.getFeatures(request)

        selectedRows = []
        for feature in selectedFeatures:
           selectedRows.append(feature.id())

        vlayer.setSelectedFeatures(selectedRows)

        output = outputBase + "/" + railway[0] + "_" + database + ".tab"

        QgsVectorFileWriter.writeAsVectorFormat(vlayer, output, encoding, coordsys, "MapInfo File", True)
        print "Translation of " + railway[0] + "_" + database + ".tab successful."


texts: Using QGIS with GDAL installed conda on Mac OSX

texts:
I used conda to install GDAL in a virtual environment in OSX 10.15 (Catalina). I have no problem running the GDAL scripts in terminal (such as gdal_translate, gdalwarp, etc).
I also installed QGIS (3.14.16-Pi). But I get an error when trying to run a QGIS command that calls GDAL.
For example, clipping a raster to a shape file using Raster > Extraction > Clip by Mask Layer yields the errors:
/bin/sh: gdalwarp: command not found
Following this post I prepended the path to gdal in QGIS: Preferences > System > Environment.  But the problem remains.
Is there some other QGIS parameter I need to set to get it find gdal?
Installing GDAL:
conda create -n gdal_env
conda activate gdal_env
conda install gdal

gdalinfo --version
GDAL 3.0.4, released 2020/01/28

Path to GDAL prepended in QGIS:
/Users/MF/miniconda3/envs/gdal_env/bin


texts: Creating Drive time analysis service area using ArcGIS Network Analyst?

texts:
Although reasonably experienced with ArcGIS I am a bit of a newbie with Network analyst.
I have a situation in which things work Ok with the tutorial data but not with mine!

I am trying to determine an 8km road distance from a set of hospitals in East Africa. I have some roads data that I can build into a network. As far as I can see there is full connectivity. If I work using detailed polygons  (which is what I want) I get this.

So the "service area does not go along roads A,B and C, which I would expect it to. and how on earth does it think I am going to get to Point D!!.

Using generalised polygons I get this:

Now I have lost the slightly anomalous point D, but the services area still does not traverse roads A and B.

I am using ArcGIS 10.2 for Desktop.


texts: Errors executing certain algorithms in QGIS Processing Toolbox?

texts:
When I execute certain algorithms in Proccesing Toolbox, for example "Convert polygon/line vertices to points" I get the following error:

  Algorithm Convert polygon/line vertices to points starting...
  shapes_polygons "Convert Polygon/Line Vertices to Points" -SHAPES "C:\Users\8725~1\AppData\Local\Temp\processing\159c5dd99017473d89a0bd321b8d8c0c\OUTPUTLAYER.shp" -POINTS "C:\Users\8725~1\AppData\Local\Temp\processing\7d7d0235a7cc4bfb993f52aa3775dfc0\POINTS.shp"
  'ascii' codec can't encode characters in position 23-33: ordinal not in range(128) See log for more details

On the other hand, when I execute the "Extract nodes" algorithm I don't get any problems...

I run on Win 10 Home 64-Bit Greek and QGIS 2.10 Pisa.
Probably the problem has to do with greek letters somewhere in my system. My user account on my pc is in greek.


texts: Projection offset between layers: OSM, GeoJSON

texts:
Don't know why there is an offset: The coordinates of the points that appear on the example were gotten through google's geocoder from some Rosario, Argentina directions. I have checked the lats and lngs and they were properly matched by the geocoder. Since I know Google works with EPSG:900913 and I've done the corresponding transformation. I cannot understand why the points are not been situated correctly.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=0">
    <meta name="apple-mobile-web-app-capable" content="yes">
    <title>OpenLayers GeoJSON Reprojected Example</title>
    <link rel="stylesheet" href="../theme/default/style.css" type="text/css">
    <link rel="stylesheet" href="style.css" type="text/css">

  OpenLayers GeoJSON-Layer as BaseLayer


texts: Help interpreting Moran's I and Geary's C results

texts:
I've got some doubts about how to interpret some results for a work. And I) would like someone could confirm or improve my conclussions.

I made a correlation analysis with SEXTANTE software in gvSIG. I obtained this results for 2.000 m and 3.000 m.

The behaviour is more or less the same in both cases, but what I want to know is how to interpret the break point at 50.000 m distance, where both index reach the higher and lower values... beyond this point, should I consider the rest of values as noise... or erratic behaviour and despise this data?

Another Question Wich distance would be nice for this type of analysis?
For this time I made an nearest neighbor analysis and I took the mean distance... Is this correct or should I take another method to decide this?


texts: Modifications to labeling engine is not supported in the QGIS composer?

texts:
I'm using QGIS 2.0 but this problem was also present in QGIS 1.8.

I have a truetype font of freeway symbols that I use to label polylines. Each symbol requires between three font characters to display correctly. For this reason, each character has a field and I have duplicated the shapefile three times and use a different field for the label each time and overlap the labels.

In the labeling engine, we can change the number of label candidates to one to force labels to overlap and this works great in the QGIS map canvas (working environment - see picture 1 below) but it's not working in the composer. The composer shows the three labels side by side instead of overlapping (see picture 2 below).

In QGIS 1.8, with labeling from the layer properties menu (that is now deprecated if I'm not mistaken) worked fine with overlapping labels but it would create "static" labels. The labeling from the layer menu creates "dynamic" labels that move around with zoom level.

I would like to have overlapping dynamic labels that also work in the composer.
How can I make the labeling engine do this?


texts: Cannot export IActiveView to PDF, invalid pixel bound value(s)

texts:
I understand that you have to re-orient the XY/min/max to fit a tagRECT but I'm still getting an "Invalid pixel bound value(s)" error. I have also tried setting Ymin/Ymax = bottom/top but I get the same error.

IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
IActiveView actView = mxdoc.ActiveView;

IExport exporter = new ExportPDFClass();
exporter.ExportFileName = @"U:\FILEPATH\topoexport.pdf";
exporter.Resolution = 300;

IEnvelope pixBounds = new EnvelopeClass();
int hdc = exporter.StartExporting();
pixBounds.XMin = actView.ExportFrame.left;
pixBounds.XMax = actView.ExportFrame.right;
pixBounds.YMin = actView.ExportFrame.top;
pixBounds.YMax = actView.ExportFrame.bottom;
exporter.PixelBounds = pixBounds;
tagRECT rectangle = actView.ExportFrame;

actView.Output(hdc, (int)exporter.Resolution, ref rectangle, null, null);
exporter.FinishExporting();
exporter.Cleanup();


texts: Updating PostGIS layer with the result of intersection

texts:
I have a PostGIS layer which I would like to intersect with another layer (adjust its borders). The result of the intersection should then become the new PostGIS-layer (update the input layer).
I could save the result as a new PostGIS table, but how could I update input PostGIS layer?


texts: Importing values into n:m relationship class of geodatabase?

texts:
I have an m:n relationship class in an ArcGIS 10.0 Geodatabase.

The values for the association table already exist in another table.

I would now like to import this table into the ArcGIS Relationship class instead of having to edit the table manually.

I don´t however see an import function.

Is there a way to import the values without manually editing the table?


texts:
I am using arcmap 10.2 . I need to find out the depth at the node points of my "point" shapefile. For this, i have tried using extract multivalues from points(Spatial Analyst).I had a .csv Events file as input point feature and bathymetry file as raster feature but it fails with an error message:

  output point feature "empty"

Does anyone know how to solve this issue?