texts: Generate point to point routes (Fiber Network not FTTH/P)

texts:
I am new to QGIS. I have a pre-existing network which I have imported into QGIS showing lines (Fiber) and points (Nodes). I am trying to put a system in place to automatically generate routes using my existing network.

At the moment I just want to generate point to point routes using my Network as the roadmap and to set up the conditions on how it routes between sites, considering each cable goes to a specific node.

Example. Cable 1 going from Node A to B (Shown in picture), then cable 2 going from Node B to c. If I trace this cable it takes the shortest path possible (Node A to C). What I am trying to achieve is when I trace this cable that it will follow Node A to Node B then to Node C.


texts: Where do 'line and point' holes come from?

texts:
It's quite common to find these topology errors when dissolving data:

Census tracts dissolved by neighborhood using ST_Union (PostGIS)

It's quite simple to understand the concept of a polygon donut hole in these cases: when the geometry really had a gap. But when your resulting holes are lines or points, my logical thinking get stuck.

What exactly makes this happen or not?
If the polygons share the same boundary lines without gaps, how does it create a line hole?
A point feature within this 'line hole' wouldn't intersect with the polygon feature?

I'm looking for a more detailed explanation, rather than the general and common "This happens when working with bad quality data"


texts: Maintaining CAD drawings and Esri Geodatabase simultaneously?

texts:
I'm wondering if there is a best practice for maintaining both CAD and ESRI GDB sytems simultaneously?

I work for my local municipality and as it is now the Eng Techs update reference drawings in CAD then I go through the changes and update the GDB. There is a lot of duplication and occasionally things get missed or fall to the wayside.

How does your organization manage 2 separate systems? Would it be better to train the Eng Techs in Esri or is it possible to streamline the data management?


texts: Using gdalwarp with multiple elevation models

texts:
I'm trying to use gdal_warp to slice and convert some SRTM data. The data I'm interested in are in S30E152.hgt and S30E153.hgt. What I'd like to do is:

Convert these DEM files from EPSG:4326 to EPSG:3857
Select a subsection of the data at coordinates 17008000,-3590694  17051000,-3503623
Write this subsection to a GeoTIFF file for further processing

Firstly, I add the two DEM files to a VRT:

gdalbuildvrt combined.vrt S30E152.hgt S30E153.hgt

Then I check the VRT to ensure the height values are sane:

gdalinfo -mm combined.vrt
[snip]
Band 1 Block=128x128 Type=Int16, ColorInterp=Undefined
    Computed Min/Max=-14.000,1501.000
[snip]

Then I use gdal_warp to convert and slice the VRT:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -te 17008000 -3590694 17051000 -3503623 -r bilinear combined.vrt sliced.tiff

Then I check the GeoTIFF to ensure the height values are still sane:

gdalinfo -mm sliced.tiff
[snip]
Band 1 Block=449x9 Type=Int16, ColorInterp=Gray
    Computed Min/Max=0.000,0.000
[snip]

It appears that whenever I use gdalwarp with a VRT containing multiple DEMs, it always clobbers the height values. For comparison, I created a VRT with only a single DEM and tried the same operations (using a smaller subsection that fits within range).

gdalbuildvrt test.vrt S30E152.hgt
gdalinfo -mm test.vrt  (produces Min/Max=-14.000,1501.000)
gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 -te 16931600 -3388300 16976200 -3439400 -r bilinear test.vrt test.tiff
gdalinfo -mm test.tiff
[snip]
Band 1 Block=446x9 Type=Int16, ColorInterp=Gray
    Computed Min/Max=53.000,1471.000
[snip]

This seems to work OK. Maybe the problem was with my coordinates. So I thought I'd warp the entire VRT, and then select the subsection:

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:3857 combined.vrt converted.tiff
gdalinfo -mm converted.tiff
[snip]
Band 1 Block=2328x1 Type=Int16, ColorInterp=Gray
    Computed Min/Max=-14.000,1501.000
[snip]
gdaltranslate -projwin 17008000 -3503623 17051000 -3590694 converted.tiff sliced.tiff

Although I've double-checked the coordinates, I always get this warning:

Warning: Computed -srcwin 914 1337 449 910 falls partially outside raster extent. Going on however.

Now this gives me a file with reasonable height values:

gdalinfo -mm sliced.tiff
[snip]
Band 1 Block=449x9 Type=Int16, ColorInterp=Gray
    Computed Min/Max=0.000,369.000
[snip]

However when I check to see if the data is sane, I get an empty image:

gdal_translate -scale 0 369 0 255 -outsize 100 200 -of PNG sliced.tiff sliced.png

At this point, I'm totally lost. I think I have data, but am not sure.


texts: Developing intersection widget using script ArcPy and geoprocessing service?

texts:
Is it possible to develop an intersection widget with geoprocessing service using ArcPy script?

I try to do that with ModelBuilder but after publishing, in geoprocessing service i don't find my inputs among the parameters.


texts: Set visible/ratio for make feature layer based on a list

texts:
So, I have a feature class with three integer fields. I need to create a feature layer based on this feature class to use it later for running some GP tools where the values of the integer fields will be used for some calculations (with ratio).

fieldA fieldB fieldC
  25     60     70

However, I don't know what field should be used, this is something that I read from a look-up table. So, when I have read the look-up table, I get a list that looks like:

 [fieldA,fieldB]

Great, now I know which fields I want to use ratio for, i.e., which one will be taken into analysis. However, I am not sure how I should approach this, that is, how can I specify in the fieldinfo variable what field should have VISIBLE RATIO based on the input list with their names. The whole fieldinfo variable is just a string and I can't see any smart way of "getting VISIBLE RATIO right after the fieldA and the fieldB" (as in the example scenario).

fieldinfo = "fieldA fieldA VISIBLE RATIO; \
             fieldB fieldB VISIBLE RATIO; \
             fieldC fieldC HIDDEN NONE"

arcpy.MakeFeatureLayer_management(input_fc,inputfeaturelayer,"#","#",fieldinfo)

I looked at the fieldinfo arcpy class, but not really sure if this will help me here and if yes how I should use it.

EDIT: before using the field_info I was forced to use the str.replace() to update the fieldinfo object.

for field in fields_list:
        updated_makefeaturelayer_fieldmapping = updated_makefeaturelayer_fieldmapping.replace('{0} {0} HIDDEN NONE'.format(field),'{0} {0} VISIBLE RATIO'.format(field))


texts: Creating centerline from lines using QGIS

texts:
How do I create a centerline (vector data) from two parallel polylines (vector data) in QGIS?
For example, if the distances between two parallel lines are 14m, then I want to create a centerline at 7m from both lines.


texts: Text color in the legend of the print composer

texts: Can I use different colors on the legend text for the title, group and object? If yes, how do you do it? My QGIS version is 2.14.12-Essen.


texts: GeoExplorer hyperlinks not working within attribute table

texts:
I'm testing the OpenGeo suite as a web mapping solution.

I have hyperlinks within the attribute table of some of my feature layers that should open up an image on my server. In my pop-up window it just shows as plan text(I can't click on it), as well as in the table itself within the query window of GeoExplorer.

I'm positive the link path itself is good. I feel like I'm missing something specific to GeoExplorer/GeoServer when dealing with hyperlinks in feature attribute tables.

Is special formatting required or maybe some settings in GeoServer?


texts: Get current route segment from Route type using ArcGIS Runtime LocationListener?

texts: I am trying to get the route segment on which I am currently navigating. i.e. I want to get the current route segment when onLocationChanged method of my LocationListener is called. A route is constituted from small segments. And we have the routing directions related to a segment but I do not know how to get the current segment in my LocationListener.


texts: Lost Foreign key during conversion of SQLite to PostgreSQL using ogr2ogr

texts:
I want to centralise several small sqlite databases into a single, big, postgresql database.  I am using ogr2ogr to convert my sqlite / spatialite (v3.7.17 / v4.1.1) databases, each one consisting of one geometry table and several attribute tables (linked to the geometry table via foreign keys), to the postgresql / postgis database ( v9.3 / v2.1.3).  Here is the code:

ogr2ogr -append -lco LAUNDER=no -lco SCHEMA=public -lco GEOMETRY_NAME=geom -f PostgreSQL PG:"..." -a_srs EPSG:4269 c:\Mydata1.sqlite line attrib1 attrib2

All my int4 primary keys (ex:  line_pk, attrib1_pk, attrib2_pk) are converted to serial type and their names are all changed to the same name: ogc_fid (eg I have ogc_fid as PK in my line, attrib1 and attrib2 tables).  I can live with this, but I find that all my foreign key constraints (from the sqlite database) are eliminated during the conversion process!  Any idea why this is so and how to correct this?  I can always re-code the constraints manually, but then I will have to somehow manage what value is put into my foreign keys each time I append a new set of data...


texts:
I would like to display all features in a feature.

Using GetFeatureInfo how would I write code to display small dots on a WMS layer, which displays pop-up and info when you click on it?

In my attribute table there is only a geometry column, there is no lat and long column. So is it possible to use GetFeatureInfo with this caveat?

I have written the following code:

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

I got this https://astuntechnology.github.io/osgis-ol3-leaflet/leaflet/05-WMS-INFO.html link for GetFeatureInfo from the Internet, but am unsure how to use this code in my code properly.

I tried to integrate it, but got a blank window as a result.


texts: openlayers 4.2 400 Error Code (Bad request)

texts:
I'm trying to load a WMS in my application, but I get the following error:

  400 (Bad Request)

Confusingly enough I still get an empty image.

The corresponding request is:

http://sg.geodatenzentrum.de/wms_dop40?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image/png&TRANSPARENT=true&LAYERS=rgb&WIDTH=256&HEIGHT=256&SRS=EPSG:4326&STYLES=&BBOX=112.5,33.75,123.75,45

I've tried setting layer projection parameters to 3857, 4326, 900913 and 25832 but none worked. QGIS also doesn't seem to like the WMS either.
So I tried simple getMap requests in my browser and the following worked:

http://sg.geodatenzentrum.de/wms_dop40?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetMap&FORMAT=image/png&TRANSPARENT=true&LAYERS=rgb&WIDTH=256&HEIGHT=256&SRS=EPSG:4326&STYLES=&BBOX=5.625,50.625,11.25,56.25

The only difference I could find was a different BBOX parameter.

Does Openlayers create a wrong BBOX parameter or what do I have to change in my code?

var view = new ol.View({
    center: ol.proj.fromLonLat([10, 51]),
    zoom: 6,
    maxZoom: 11,
    projection: 'EPSG:3857'
});

var luftbilder = new ol.layer.Tile({
    source: new ol.source.TileWMS({
        url: 'http://sg.geodatenzentrum.de/wms_dop40?',
        params: {'LAYERS': 'rgb',
                'VERSION': '1.1.1',
                },
        projection: 'EPSG:4326'

    })
});

var map = new ol.Map({
    target: 'map',
    view: view,
    layers: [luftbilder]
});

weird, fiddle works fine fiddle (in Chrome)


texts: How to write XML attributes from feature with ArcView license?

texts:
I need to generate xml of attribute data of a feature, but I have only Arcview license so, I can't use Export XML tool. What are my options?

What is wrong with this code (after applying search cursor writing the attributes of node):

 for row in rows:
  order.attrib['name'] = row.getValue('Id')
  order.attrib['UTM'] = row.getValue('UTM')

I am not able to see any addition to order tag. missing something...
Thanks

Ok - just realized there was some issue with the feature:
Here is the code and error message after the change (So something is broken fundamentally):

>>> import xml.etree.ElementTree as xml
>>> import arcpy
>>> fclass= arcpy.SearchCursor(r"..\SiteMarker_Project")
>>> order = xml.Element('Order')
>>> site = xml.Element('Site')
>>> order.append(site)
>>> rows = fclass
>>> for row in rows:
  order.attrib['name'] = row.getValue('Id')
  order.attrib['UTM'] = row.getValue('UTM')

Traceback (most recent call last):
  File "<pyshell#8>", line 3, in <module>
    order.attrib['UTM'] = row.getValue('UTM')
  File "C:\Program Files\ArcGIS\Desktop10.1\arcpy\arcpy\arcobjects\arcobjects.py", line 1011, in getValue
    **return convertArcObjectToPythonObject(self._arc_object.GetValue(*gp_fixargs(args)))
RuntimeError: ERROR 999999: Error executing function.**


texts:
I have a shapefile to which I joined two tables from a MySql database. I used one column of the first table to color the polygons in my shapefile which works fine. Three columns out of the second table should now be used to produce a pie chart on top. I can easily select those features, so the join worked. However, no pie chart is drawn. "Display diagrams" is checked.
I tried the following so far:

Change diagram type (Did not work)
Try to display a value from the first table instead of the desired table. (Did not work)
Changing the priority from low to high. (Did not work)
Delete the first table.(Did not work)
QGIS is not drawing anything.

The log for errors and warnings are completely empty. Does anyone has an idea? Thanks in advance.


texts: Alternative indexing methods for point set operations

texts:
It is common to use a bounding box spatial index to improve performance when working with a large number of features.  Where operations are performed against individual geometries with a large number of vertices do any similar optimization strategies exist?

For example, do any data structures exist that can speed up point in polygon, or union operations?


texts:
I try to filter multiple values in the same column in CartoDB. I then like to filter on another column.

I do not manage to filter on multiple values in one column in the graphical interface (one value ok, but more I see no way to do this).

I then write the SQL query:

    SELECT * FROM bdi_adm2
WHERE name_1 ILIKE '%bururi%'
OR name_1 ILIKE '%cancuzo%'
OR name_1 ILIKE '%karuzi%'
OR name_1 ILIKE '%makamba%'
OR name_1 ILIKE '%rutana%'
OR name_1 ILIKE '%ruyigi%'

But then I am not able to filter on another column (it tells me that I have to clear the selection first).

How can I filter through the input box for multiple values?


texts: Split polylines at polygon boundaries in QGIS

texts:
I have a set of contour lines and a polygon set containing glaciers.
Because I want the glaciers to have blue contours, I need to split the contour lines at the glacier boundaries.

I can achieve this using the following method:

1) Intersect the contours with the glacier polygon to create the glacier contours.
2) Create a spatial difference between the glacier polygon and an extent polygon.
3) Intersect the contours with the spatial difference from (2) to create the blue conoturs.
4) Merge the two results together if I want it in a single data set.

But the following method would be more convenient, and usable for other types of analyses as well.

1) Combine a polyline dataset (or point dataset for that matter) with a polygon dataset.  The lines with be split at the polygon boundaries and each output line (or point) will take the attributes from the polygon it is inside.

Such a function would make my problem so much simpler to solve, as I can do it in one single operation.

Is there a function for this in QGIS?


texts: How to enter custom string values into GPString multiValue parameter in 10.4 python toolbox?

texts:
I am trying to be able to add a list of years (or return periods) to a python toolbox tool I want to create. In the getParametersInfo function, I am currently setting up a parameter as a "GPString" datatype and want it to be a multiValue list. Here is my current code:

import arcpy
class Toolbox(object):
    def __init__(self):
        self.label = "Toolbox"
        self.alias = ""
        # List of tool classes associated with this toolbox
        self.tools = [TestParameters]
class TestParameters(object):
    def __init__(self):
        self.label = "TestParameters"
        self.description = ""
        self.canRunInBackground = False
    def getParameterInfo(self):
        param0 = arcpy.Parameter(
                    displayName ='Do Nothing Here:',
                    name ='in_features',
                    datatype ="GPFeatureLayer",
                    parameterType ='Required',
                    direction ='Input')
        param1 = arcpy.Parameter(
                    displayName ='Enter Return Periods as Years',
                    name ='tbl_features',
                    datatype = "GPString", #"GPValueTable", #"GPTableView", #"DETable", #
                    parameterType ='Required',
                    direction ='Input',
                    multiValue=True)
        param1.filter.type = "ValueList"
        param1.filter.list = []
        params = [param0,param1]
        return params
    def isLicensed(self):
        return True
    def updateParameters(self, params):
        param1.filter.list = []
        return
    def updateMessages(self, params):
        return
    def execute(self, params, messages):
        arcpy.AddMessage("Testing Tool, nothing here yet")
        arcpy.AddMessage(params[1].valueAsText)
        return

However, when I pull up the tool in ArcCatalog and open the tool prompt I get an empty box. I expected that, but the "Add Value" button does not allow me to add any values. I would like to add years values such as "2", "5", "10", "100", "500", & etc...

Am I missing something that should go inside the updateParameters function? Or is there a way to enable the "Add Value" button? It would be nice not to have to hard code year values, as others may be using this script later and may have different years they need to input. This is why I left the param1.filter.list = [] as a blank list, as the ESRI documentation for it doesn't cover how to enter values after the tool is opened.

Or is there a different type of parameter that can easily do this? I've tried the values commented out on the datatype line for param1. If any of those would work, which one might help me be able to customize this?


texts: Is there a reason why my ArcGIS Server token works, when the URL refers to the physical server, but does not work when I use the mapped URL?

texts:
http://servername:6080/arcgis/rest/services/xxx/xxx/MapServer?token=xxx

works fine and I am directed to the service information without a need to login

however

when using the mapped (external DNS record) reference, it gets to the ArcGIS REST login page

http://xxx.com.au/arcgis/rest/services/xxx/xxx/MapServer?token=xxx


texts:
I have created a topology from OSM data and would like to simulate a car driving every possible route along the map. The purpose is to "survey" the osm data as if in a car. My ultimate objective is to find every unique journey segment of five miles, I am aware that there will be a LOT of them.

I am using pgrouting/postgis and have created a topology from roads. I've also written an algorithm that given any osm_id  will find all connected osm_id's to a depth of N. This works by finding all connected osm_ids and then recursing through these until the depth is reached.

I could modify this to change N to a distance of five miles.

But is there a better way of traversing every possible path in a network?

For information, here is my current algorithm which is sub-optimal as it does not deal well with junctions that have more than one connection

** Find connections **

CREATE OR REPLACE FUNCTION ng_research.get_connections(id bigint, table_param character varying)
RETURNS bigint[] AS
$BODY$
DECLARE
  ret_vals integer [];
BEGIN
--quite simply given a source node for a line string returns all of the connected nodes including its own target

WITH CTE AS (
 SELECT source from ng_research.all_roads WHERE osm_id=id
)

SELECT array_agg(osm_id)  into ret_vals FROM ng_research.all_roads AR,CTE  WHERE  AR.target=CTE.source or AR.source=CTE.source;
RETURN ret_vals;
END;
$BODY$
 LANGUAGE plpgsql VOLATILE

** Build linestring arrays **

CREATE OR REPLACE FUNCTION ng_research.road_string_finder(node bigint, table_param character varying, depth integer default 0, max_depth integer default 2, previous_nodes bigint[] DEFAULT NULL::integer[]) RETURNS bigint[] AS
$BODY$
DECLARE
 connected_nodes bigint[];
 node_var integer;
BEGIN

IF depth = 0 THEN
 DROP TABLE IF EXISTS ng_research.nodes_visited;
 CREATE TABLE ng_research.nodes_visited (id bigint);
END IF;

previous_nodes := array_append(previous_nodes, node);

SELECT ng_research.get_connections(node, table_param) INTO connected_nodes;

IF connected_nodes IS NOT NULL THEN
 FOREACH node_var IN ARRAY connected_nodes LOOP

  IF NOT (node_var = ANY(previous_nodes)) THEN
   RAISE NOTICE '% :  %', node_var, previous_nodes;
   IF depth < max_depth THEN

    IF (SELECT node_var from ng_research.nodes_visited WHERE id=node) IS NULL THEN
    SELECT ng_research.road_string_finder(node_var, table_param, depth + 1, max_depth, previous_nodes) INTO previous_nodes;
    INSERT  INTO ng_research.nodes_visited (id) VALUES (node_var);
   END IF;
  END IF;
 ELSE
 RAISE NOTICE 'SEEN YOU';
END IF;
END LOOP;

END IF;

IF array_length(previous_nodes,1) > 1 THEN

 previous_nodes := array_remove(previous_nodes, node);
 END IF;
 return previous_nodes;
 END;
$BODY$
LANGUAGE PLPGSQL;

** find connections for a single linestring **

SELECT ng_research.road_string_finder(204359430, 'ng_research.all_roads',0,5);


texts: Is it better to take the elevation of a GPS unit or that of a Digital Elevation Model?

texts:
When searching of gps routes in wikiloc, for the trailrunning race "Transvulcania" I have several routes to look. I downloaded a lot, but every route has a distance between 71 kilometers and 74 kilometers, of course someone in the race could be lost for a while, so that is not a big problem for now, but when I see the elevation gain, is very different between them.
For example the route in wikiloc says:
74.460m +4.171m -3.874m
The same route in eMapzone
71.923m +6.495m -6.132m

Another one
wikiloc
72.100m +4.560m -4.210m
eMapzone
72.045m +5.638m -5.286m

and the worst is that the official distance of the race is: 73.000m and +4.300m.

So I decide to normalize the elevation of the route with a Digital Elevation Model, I upload the original file to http://www.gpsvisualizer.com/elevation and then I look the information with eMapzone, and there is some consistency in the data.

The elevation DEM calculated of the same route with several and diferents gps hardware is something in between +5354 , +5598 and -5015, 5263.

So my question is? is better to take the elevation of a Digital Elevation Model?


texts: Assigning irregular points to a grid and assigning grid number as attribute

texts:
I have a large point data set that has irregular geometry (random sampling) that I want to bin in 2D XY space. I've attached a sketch below to illustrate what I mean. Apologies for hand drawing, hope no one minds.

So far I've been using the fishnet option in ArcGIS+QGIS and SAGA tool polygon attribute to points. However I have reached the workable file size limit for my shapefiles (>2GB). When attempting to use feature classes ArcGIS is incredibly slow at performing the required task, so I have been mostly focusing on SAGA for its speed.

It's been suggested to me by a colleague that I could use GMT's blockmean, blockmedian, blockmode programs. However I have limited experience in this but would be willing to learn IF it is the right direction

I can see the maths of all this as a matrix operation, but the problem I am having is with the large data size. My point files are >700mb in just XY coordinates (billions of points).

I am not looking for interpolation, just looking to assign the points a 'bin' reference in XY space.


texts: arcobjects 10.3.1 c#: Loading Addins programmatically

texts:
I am writing addins in C#, and am trying to load an addin. I thought i had the solution, but it seems like the full code is not loaded.

The button is added to my toolbar, and code seems to run, too a point. Once i need to use the singleton extension of the child, it crashes.

The child addin will run by itself, with no issues.

Calling Addin (Parent):

public class ParentExten : Extension
{
 .....

/// <summary>
        ///
        /// </summary>
        protected override void OnStartup()
        {
            WireDocumentEvents();
            // load other extensions:

            // Call FindExtension to load this just-in-time extension.
            var ais =  ArcMap.Application.FindExtensionByName("ChildExtension");

        }

}

Child Addin: config.esriaddinx

<AddIn language="CLR4.5" library="name.dll" namespace="name">
    <ArcMap>
      <Extensions>
        <Extension id="childExtension" class="classname"  showInExtensionDialog="true" autoLoad="false" />
      </Extensions>
      <Commands >
        <Button  id="button_id" class="namspace.classname"
                  image="Images\Layer.jpeg"  onDemand="true"/>

        <Tool id="tool_id" class="hiddenname"
             category="Layer" image="Images\Layer.jpeg" onDemand="true" />
      </Commands>
    </ArcMap>
  </AddIn>


texts: Shortest path analysis has unattainable vertices in result from PyQGIS

texts:
You can try this code directly in the console. It basically finds the shortest path but I only have unattainable vertices in the result.

Do you know how I can correct this?

My layer consists of MultiLine geometries.

layers = iface.mapCanvas().layers()
segmentLayer = layers[0]
loggerLayer = layers[1]

features = segmentLayer.getFeatures()
# self.features(features)

director = QgsVectorLayerDirector(segmentLayer, -1, '', '', '', 3)

crs = segmentLayer.crs()

transform = QgsCoordinateTransform(loggerLayer.crs(),crs, QgsProject.instance())

builder = QgsGraphBuilder(crs)

# graph = builder.graph()
pStart = QgsPointXY(677725.236666,6155028.577216)
#pEnd = QgsPointXY(676933.041661,6154972.539126)
director.addStrategy(QgsNetworkSpeedStrategy(10,1.0,1.0))
tiedPoints = director.makeGraph(builder, [pStart])
graph = builder.graph()

startId = graph.findVertex(tiedPoints[0])
#endId = graph.findVertex(tiedPoints[1])

(tree, cost) = QgsGraphAnalyzer.dijkstra(graph, startId, 0)

P.S. The shortest path algorithm is working in the toolbox (predefine network analysis functions)


texts: Comparing type of fields of different layers in Qgis

texts:
Is it possible to have two property layer dialogs open at the same time? I would like to see the properties of two layer's fields at the same time.
I need this because when I want to join two layers by a common field I like to check they are of the same type. I am using 2.18.11.


texts: v.generalize in QGIS throwing "Unable to create table" error

texts:
I am trying to smoothen few polyline shapefiles using v.generalize in QGIS 3.4. But it's not working. Even i tried the same in QGIS 2.18. Getting the same error.

Tried the same process in GRASS 7.4.2. The shapefile is not even getting loaded. When v.in.ogr is used, getting the exact "Unable to create table" error


texts: ArcGIS10.2 SDE relationship class referencing a deleted feature - how to clean up relationship class?

texts:
I am looking for a smart approach to cleaning up a relationship class, which consists of relationships that look at previously deleted objects.
To illustrate, we have db.schema.SITE (spatial), db.schema.STATUS (non-spatial) tables in our multi versioned SDE geodatabase (SQL Server), and a relationship class that relates the two tables called db.schema.SITE__STATUS. However, when I run the following query, it returns over 4000 records:
select
ss.*, status.*, site.*
from
db.schema.SITE__STATUS ss
left join db.schema.SITE site  on site.OBJECTID = ss.SITEID
left join db.schema.STATUS status on ss.STATUSID = status.OBJECTID
where
ss.statusid not in (select objectid from db.schema.STATUS)
and
ss.hubsiteid  in (select objectid from db.schema.SITE)

I suspect over time statuses have been deleted but somehow the relationship class failed to get updated. But I don't understand (and would love to understand) how this situation is allowed to happen in the geodatabase in the first place.
One approach to cleaning up would be to remove the returned records from the relationship class using SQL, but I fear this will be filled with GOTCHAs I failed to foresee.
What would be a smart way to approach this cleaning up task?


texts: Creating grid from spatio temporal point data?

texts:
I have a spatio-temporal point data of species sighting with monthly observation for random locations within my study area. I wish to create a grid(fishnet) for each month for each year. Also, certain grids for a month and year has multiple points so need to derive mean.

Old [I have point data for an area over a period of years (2010-2014). I would like to create a grid of 0.25 X 0.25 decimal degrees for each month. Also to derive a mean where there are more than one point in a grid mesh for that particular month. Purpose of creating this grid is to extract the values from rasters based on the grid.]


texts: QGIS Layer Order Panel - add layers at top of layer order

texts: I have a QGIS workspace in which the order of the layers is controlled by the Layer Order Panel.  When I add a layer it automatically ends up at the bottom of the order and has to be manually shifted to the top so that it's visible.  Is it possible in the options to specify that added layers end up at the top of the layer order instead of the bottom?


texts: Filter a collection using nightime in Google Earth Engine

texts:
I'm working with GOES-16 full disk data, trying to calculate the frequency of low clouds in northern Chile.
I need to filter the collection so it contains only images from September 2017 at night. This means, between 22 UTC and 11 UTC. To filter the year and month I have used ee.Filter.calendarRange but I don't know if it is working with the night hours because the start argument (22) is greater than the end (11).
Here is the code I have been using:
//creation of temporal filters
var rYear = ee.Filter.calendarRange(2017,2017,'year'); //year 2017
var rMonth = ee.Filter.calendarRange(9,9,'month'); //september
var rHour = ee.Filter.calendarRange(22,11,'hour'); // nocturnal hour range in UTC (local time: UTC-3)

//filter collection by temporal range
var night = GOES16.filter(rYear).filter(rMonth).filter(rHour);

Here is what the reference page of ee.Filter.calendarRange says, but I don't understand it:

If end < start then this tests for value >= start OR value <= end, to allow for wrapping.

And here is the complete code:
https://code.earthengine.google.com/27d18828f97161552d812cf82aa87477?noload=true
Is there a way to filter the collection only in the desired hour range (22-11 UTC) for a whole month and for a specific year?


texts: Finding distance and time from multiple points to nearest destination(s) using ArcGIS Network Analyst?

texts:
I am trying to compute the distance and time taken for each origin (40,000+ house points) to the nearest train station (there are 111 of them). How can I do that in ArcGIS?

So example Houses H1 to H10 nearest station is S1, and H11 to H30 is S2, so on...

Update: The distance would be network distance not straight line.


texts: Geometry geografic transformation with OGR, OSR - UTM to Degree

texts:
I am trying make a geotransformation from UTM coordinates to decimal degree coordinates, but i have only zeros as result... I must be missing something, i already take a lock on other users codes but i can't find my bug. I would apreciate a little help here!
Those are the shapefiles i used (i uploaded in Google drive):
https://drive.google.com/drive/folders/1URP8UxjNxaI8XO-qNrzWb8-lNA1ZDCBe?usp=sharing
sentinel_grid_shp_path = r'grid_s2.shp'
shp_path = r'zonal_3_zones_1000_m2_filtering_L4.shp'

from osgeo import ogr, osr

# OSGEO Driver
driver = ogr.GetDriverByName("ESRI Shapefile")

##### datasource 01 - 4674
dataSource = driver.Open(sentinel_grid_shp_path, 0)
dataSourceLayer = dataSource.GetLayer()

# Get DATUM information
sourceDatumRef = dataSourceLayer.GetSpatialRef()
# Get EPGS code
a = sourceDatumRef.ExportToWkt()[sourceDatumRef.ExportToWkt().rfind("AUTHORITY") + 9:].replace("]]", "]")
b = a.replace("[", "").replace('"', '').replace("]", "")
### EPGS code
sourceDatumRefEPGS = int(list(map(str, b.split(",")))[1]) ### 4674

# Spatial ref 4674
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(sourceDatumRefEPGS)

##### datasource 02 - 32721
dataReceiver = driver.Open(shp_PATH, 1)
dataReceiverLayer = dataReceiver.GetLayer()

# Get DATUM information
receiverDatumRef = dataReceiverLayer.GetSpatialRef()
# Get EPGS code
a = receiverDatumRef.ExportToWkt()[receiverDatumRef.ExportToWkt().rfind("AUTHORITY") + 9:].replace("]]", "]")
b = a.replace("[", "").replace('"', '').replace("]", "")
### EPGS code
receiverDatumRefEPGS = int(list(map(str, b.split(",")))[1]) ### 32721

# Spatial Ref 32721
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(receiverDatumRefEPGS)

# Definition of geotrnaformation
geo_transform = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

### Parsing througout layer features and applying geotransformation
for feature in dataReceiverLayer:
    geom = feature.GetGeometryRef()
    centroid_wkt = geom.Centroid().ExportToWkt()
    centroid_geometry = ogr.CreateGeometryFromWkt(centroid_wkt)

    # GEOTRANSFORMATION - ONLY '0' AS RESULT!!!
    print(centroid_geometry.Transform(geo_transform))

del feature, geom, centroid_wkt

# Excuding and deleting all objects

driver = None

dataSource = None
dataReceiver = None

dataSourceLayer = None
dataReceiverLayer = None

sourceDatumRef = None
receiverDatumRef = None

sourceDatumRefEPGS = None
receiverDatumRefEPGS = None

inSpatialRef = None
outSpatialRef = None

a = None
b = None

centroid_geometry = None

ogr = None


texts:
I defined some virtual layers in my QGIS project. Running a compiled QField in Ubuntu environment, virtual layers works fine. But using QField in Android (last official from playstore or the compiled apk from Ubuntu) virtual layers in the project cannot be loaded and 'invalid data provider virtual' message appears.

Is it possible to use virtual layers with QField in Android?


texts: How to use the GeoServer REST api to get a PNG representation of a layer?

texts:
Does anybody know if I can retrieve by REST service in GeoServer an raster representation of layer?
Just like in for ex. ArcGIS REST Api:
When I send:

MapServer/export?dpi=96&transparent=true&format=png8&bbox=338477.336%2C360036.34849999996%2C619993.648%2C548154.8615&bboxSR=2180&imageSR=2180&size=1064%2C711&f=image

In response, I get png image of layer in defined bbox.


texts: Many of the Esri samples fail the WC3 HTML validator. Is this a problem?

texts:
I'm running into problems with some ArcGIS Server JavaScript code (page elements are positioned different in IE than FF/Chrome), so I checked it against the HTML validators at http://validator.w3.org and http://html5.validator.nu

Many errors were found, so I returned to the official Esri JS samples to ensure I was using valid code. Surprisingly, I found that many of the samples also failed the validators.

Examples of errors include:

Error: Attribute label not allowed on element button at this point.
<button id="dropdownButton" iconClass="bingIcon" label="Base" data-dojo-type="dijit.form.DropDownButton">

from the Basemap sample.

Is it a problem that the official Esri samples fail the validators? Is this related to the use of DOJO, and therefore not an actual problem?


texts: ArcObjects, "License not intialized" error, but it should be initialized

texts:
In my VB.NET application, I am using ArcObjects 10.4 and I am initializing ArcGIS. I load a feature class from a shapefile. The problem is that I get an HRESULT error 0x80040228, which is 2147220952 in decimal and it means "License not intialized" as you can see on this website: http://resources.arcgis.com/en/help/arcobjects-net/componenthelp/index.html#//0025000000mn000000

The strange thing is, I load two feature classes. One from a shapefile, then another from a file geodatabase. As I said, loading from the shapefile crashes, but not from the FGDB. Also, if I invert the two LoadFeatureLayer lines, the FGDB feature class loads, and then, the shapefile feature class loads too. I tried with different features classes, the result is the same.

Is it a normal behavior?

Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.Geodatabase
Imports ESRI.ArcGIS.esriSystem

Public Class frm_Loading
Public aoInit As IAoInitialize

Private Sub frm_Loading_Load(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles MyBase.Load
    InitializeArcGIS()
    LoadFeatureLayer("export_output_2.shp", "c:\temp")
    LoadFeatureLayer("C_R137_CH", "c:\temp\FGDB.gdb", True)
End Sub

Public Sub InitializeArcGIS()
    If aoInit Is Nothing Then
        If ESRI.ArcGIS.RuntimeManager.ActiveRuntime Is Nothing Then
            ESRI.ArcGIS.RuntimeManager.BindLicense(ESRI.ArcGIS.ProductCode.Desktop)
        End If

        aoInit = New AoInitialize
        aoInit.Initialize(esriLicenseProductCode.esriLicenseProductCodeBasic)
    End If
End Sub

Public Function LoadFeatureLayer(ByVal layer As String, ByVal cheminWorkspace As String, Optional FGDB As Boolean = False) As IFeatureLayer
    Dim aWorkspaceFactory As IWorkspaceFactory = Nothing
    Dim aFeatureWorkspace As IFeatureWorkspace = Nothing
    Dim aFeatureLayer As IFeatureLayer
    Try
        If FGDB Then
            Dim t As Type = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory")
            Dim obj As System.Object = Activator.CreateInstance(t)
            aWorkspaceFactory = obj
        Else
            Dim t As Type = Type.GetTypeFromProgID("esriDataSourcesFile.ShapefileWorkspaceFactory")
            Dim obj As System.Object = Activator.CreateInstance(t)
            aWorkspaceFactory = obj
        End If
        aFeatureWorkspace = aWorkspaceFactory.OpenFromFile(cheminWorkspace, 0)
        aFeatureLayer = New FeatureLayer
        aFeatureLayer.FeatureClass = aFeatureWorkspace.OpenFeatureClass(layer)
        aFeatureLayer.Name = aFeatureLayer.FeatureClass.AliasName
        Return aFeatureLayer

    Catch ex As Exception
        Debug.Print(ex.ToString)
        Return Nothing
    Finally
        If aWorkspaceFactory IsNot Nothing Then System.Runtime.InteropServices.Marshal.FinalReleaseComObject(aWorkspaceFactory)
        If aFeatureWorkspace IsNot Nothing Then System.Runtime.InteropServices.Marshal.ReleaseComObject(aFeatureWorkspace)
    End Try
End Function
End Class


texts:
I am working in QGIS 3.4.8-Madeira on MacOS 10.14. I simply want to create a buffer on a point shapefile, but the layer is not being read in UTM but in degrees. I understand this is a common question but am having trouble with solutions I have found. I want the CRS to EPSG:32642/WGS 1984 UTM 42N (Afghanistan).

My steps to change the CRS were to:

Set project CRS (EPSG:32642/WGS 1984 UTM 42N);
Import the .csv with the desired CRS;
Save the .csv as a .shp with the desired CRS.

I have tried this set of steps several times. While the .qgz and the .prj of the shapefile all state they are in EPSG 32642, my points keep rendering at 0,0, and coordinates are still in decimal degrees rather than meters. I believe the source of this issue may be in the .csv formatting. How do I convert the decimal degrees to meters? Apologies for the question repetition.


texts: PostGIS Vector Tiles zoom affects tile production

texts:
I'm building a small tile service to present postgis data on the fly to a leaflet map.

for testing I am using the data set  ne_110m_admin_0_countries

The first problem is the extent of the tiles returned, many are error'ing in predictable patterns. For example on zoom 0 no tiles below the equator are drawn.

between zoom 0 and 3 New Zealand doesn't exist

and also at zoom 3 nothing appears below 41 degrees south, however this behaviour does not exist reflected in the norther hemisphere.

I'm using the following to serve up the tiles on a flask server:

@app.route('/tiles/test/<int:z>/<int:x>/<int:y>', methods=['GET'])
@login_required
def tiles(z=0, x=0, y=0):

    tile = get_tile(z, x, y, 'ne_110m_admin_0_countries')
    if tile != None:
        response = make_response(tile)
        response.headers['Content-Type'] = "application/octet-stream"
        return response
    else:
        return ('', 502)

def get_tile(z,x,y,TableName):
    cursor = db.connection.cursor()

    query= """WITH webmercator(envelope) AS (
              SELECT ST_TileEnvelope({},{},{})   -- z, x, y
            ),
            wgs84(envelope) AS (
              SELECT ST_Transform((SELECT envelope FROM webmercator), 4326)
            ),
            geometries(id, geom) AS (
              SELECT id, ST_Transform(geom, 3857)
              FROM  {}
              WHERE geom && (SELECT envelope FROM wgs84)   -- do not read objects outside the tile
            )
            SELECT ST_AsMVT(tile) FROM (
              SELECT id, ST_AsMVTGeom(geom, (SELECT envelope FROM webmercator))
              FROM geometries
            ) AS tile""".format(z, x, y, TableName)
    try:
        cursor.execute(query)#,(z,x,y))
        tile = bytes(cursor.fetchone()[0])
        cursor.close()
        return tile
    except:
        cursor.execute("rollback")

The time zones from natural earth also show weird issues when zooming in out. For example at zoom level 4 the longitude extent of the tiles is up to 158 degrees east/west, at zoom 3 its 135 degrees, and at zoom 2 it's restricted even more to only 90 degrees east west

Any suggestions on why the zoom might be effecting the tiles shown?

Stay safe all.


texts: QGis 2.8.1 Series of maps from different layers

texts:
I have many layers in one QGis-Projekt based on the same .shp-File (Community borders)  with different fields classified and displayed in each layer (Mean Income, Mean Age, Average Houshold Size ...).

I want to use the same template in the map composer for all layers. Is there is any way to export all layers at once with the same template so that I don't have to go through all the layers individually?


texts: Install latest Osmosis version on Ubuntu

texts:
I am using Ubuntu 14.04 and I would like to know the easiest way to install the latest (stable) version of Osmosis, which is 0.44.

I ran

sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable
sudo apt-get update
sudo apt-get install osmosis
osmosis -v

That last command returns many lines, but the second is

INFO: Osmosis Version 0.40.1

0.40 seems to have too many problems with 64-bit integers. So I want to download the latest stable, but I don't know where Osmosis is currently installed to and thus where I should copy the downloaded files to.

http://wiki.openstreetmap.org/wiki/Osmosis#Downloading


texts:
What I am trying to do is to take a specific range of time (summer) over multiple years and then sum the total amount of precipitation onto a final layer.

This is the code I am working with

var dataset = ee.ImageCollection('IDAHO_EPSCOR/GRIDMET')

var precip = dataset.select('pr');
var precipVIS = {
  min: 0.0,
  max: 10.0,
  palette: ['d8d8d8', '4addff', '5affa3', 'f2ff89', 'ff725c'],
};
Map.setCenter(-86.356, 44.686, 9);

Map.addLayer(
    precip.filterDate("2014-06-15", "2014-07-15"), precipVIS,
    '2014');
Map.addLayer(
    precip.filterDate("2015-06-15", "2015-07-15"), precipVIS,
    '2015');

Now I am looking for the last part where I sum the results of the two map layers '2014' & '2015'.


texts: Collector App dropping collected points from one specific layer after submission?

texts:
I'm using the collector app for iOS (iPad Air 2) for water utility asset data collection. For every layer or asset type I can mark a point, fill out the form, attach photos, and submit. No problem. It immediately appears on the collector map and syncs with ArcGIS online. The problem is when I do this with valves and valves specifically, all of the above appears to work correctly, except the point on the map disappears after about 2 seconds. The valve attribute field is empty, none of the points, attachments, or associated fields are being saved. Everything else is perfect. Collector isn't throwing any errors and I cannot find out what's causing it to delete valve points. All of these layers are on the same feature service, created by ESRI. The valves layer has edit permissions and displays at the same scale as everything else.

Is this a common bug? Am I doing something horribly wrong? I was working in an area with excellent cell service; I never went offline. Also, only points from this one layer were dropped; the hydrant points I collected right next to these valve points went through flawlessly. Even sitting here back in the office I can create little test points and the exact same thing happens: valves disappear, everything else saves.


texts: ArcGIS raster ASCII range not displaying as expected?

texts:
On ArcGIS, I have initially used the mosaic feature to join a number of .asc files. The mosaic output displays the same range (max and min) of these .asc files. I  now need to convert the raster mosaic file to an ascii file but for some reason the output of this conversion is displaying a different range. This is a bit strange to me  (I must say that I am  not a gis expert).

Does anyone have an idea of why this is happening?


texts:
How can I get layout widgets objects dynamically? In case of typical widgets it can be achieved by

  FlexGlobals.topLevelApplication.widgetManager.getWidget()

, however i can't find a way to get to the layout widgets like navigation or map switcher.


texts: What expression will give age in years in QGIS?

texts: I have the date at which a pipe was laid. I want to create a field giving the age of the pipe in years. How do I do this in QGIS?


texts: How to simplify streets such that seperate driving tracks are reduced?

texts:
Streets in Openstreetmap have more detail than I need. For instance I do not need seperate street tracks for each driving direction, instead I just want to show the general path of the street. Example:

OSM linestrings for a street

Simpified linestring(s) (here it is just one, but in case of sidearms would be > 1)

input:  collection of linestrings
output: collection of linestrings
simplified lines = f(osm lines)
f = ?

the output is a collection as well, since sidearms should remain available.

How do I programatically simplify such linestring collections?

Would like to use JTS, NTS, GEOS but not PostGIS if possible.


texts: smallest, lossless, GDAL ready image format

texts:
I would like to know what the smallest file format is I can use with GDAL that fulfills the following criteria:

lossless compression
(relatively) easy to read/write with GDAL Python bindings
least HDD space used
acceptable read/write performance

What have I tried so far? The best result so far seems to be a tiled GeoTiff with LZW compression:

gdal_translate -co "TILED=YES" -co "COMPRESS=LZW"

Why am I looking for this:

I am about to process time-series datasets which will take up around 6-8 TB in band sequential binary format. The processing is quite extensive but the output will only weigh in at around 30-40 GB. I know I am looking at a hit in read/write performance when using compression. An input/output speed degradation of around 300%-500% is acceptable - occupying 6-8 TB HDD space is not.

edit:

After the answer I wrote a small python script which tests the performance of different compressions and an input image of your choice and a little description - GeoTiff compression comparison


texts: Is it possible to remove the QGIS icon from dialogs windows using pyqt?

texts:
I am developing a QGIS plugin using PyQt 4.x and the QGIS API. The plugin uses several dialog windows (QDialog Class). Some of them (I don't know why some do and others don't) show the Qgis icon on the top left part of the frame. I am trying to remove all icons from all dialogs using Qt.WindowFlags.

In most of them the following works:

dialog_window.setWindowFlags(Qt.CustomizeWindowHint | Qt.WindowTitleHint)

result:

when I want to keep the closewindowHintButton:

dialog_window.setWindowFlags(Qt.CustomizeWindowHint | Qt.WindowTitleHint | Qt.WindowCloseButtonHint)

the logo does not go away and I end up with this:

I also tried to set up the setWindowIcon() method to an empty QIcon() instance, still the logo does not want to disappear.

Has anyone ever found a solution to this?


texts: Looking for free Mexican transportation datasets

texts:
I have been looking for data of Mexico. What I am looking for is the data of a Transportation in Mexico, specifically Chihuahua province. I have checked out on ArcGIS online and Google but came up nothing good. Is there one that I overlooked? I am looking for free gis data.

Thank you.


texts:
I am new to OSM and I would like you to help me with some questions I have.
I'm styling my raster map of Europe.
The data for the raster I have obtained from srtm.csi.cgiar.org but the data for Norway, Sweden does not appear and I would like to obtain them. So I have obtained them from land.copernicus.eu.
Once all the data has been obtained, I have tried to make a hybrid solution to add this information, but on the one hand the data is stepped on since there is a small common area of both in the downloaded data.
Is there any way so that the information for the common area of Norway does not appear double or is there somewhere where I can download all the data for the raster of the map of Europe that will quickly load the map and thus not have to do a hybrid solution with information from two sites? I've tried downloading everything from Copernicus and it slows down and doesn't even load the map.


texts: Can't copy field to another: 'Empty value for Objectid = 0'

texts:
I tried to copy a field with dates that had type 'TEXT' to a new column of type 'DATE' because it is not possible to change the initial type of the original field.

arcpy.CalculateField_management(r'C:\Users\USER\Desktop\Deuv.shp','DC_DATE','DC_DATE1')

The initial field has all rows with data but it gives this error:

  'Empty value for Objectid = 0'

And it can't copy anything to the other field.


texts:
I followed instructions in this answer to Creating Directional Flow Arrows for sewer lines? to add directional arrows to the one-way streets on a map I'm creating.

I've currently created a marker line with the standard arrow marker but there are a lot of one-way streets (and one-way street segments). So even with just a marker on the centerpoint of each line, I get something like the below image

I'm worried that if I dissolve streets together by name, that the direction of their drawing might get muddled (and the labels point in the wrong direction).

Do tutorials exist for labelling one-way streets in QGIS?


texts: How can I extract each value of a column from a Feature Collection and append it to a new list in Google Earth Engine?

texts:
I want to map a function over a Feature Collection to extract some column values and turn it into a list.
I do know how to map a function over collections, but I'm new with Google Earth Engine and I don't know very well javascript, so I don't know how to work with lists and other features from this language.


texts: Is it possible to have a legend and other 'decorations' above the raster image in geoserver?

texts:
I would like to use geoserver to subsample a global raster geotiff into many smaller, regional images. These images all need to have their own legend and descriptive text. In all of the examples that I have seen with geoserver, the legend and text overlay and partially obscure the actual raster image.
Examples: http://docs.geoserver.org/latest/en/user/services/wms/decoration.html#wms-decoration-types, http://geoserver.geo-solutions.it/edu/en/pretty_maps/map_decorations.html

What I would like is a result like this:

I know that technically the legend graphic can be requested separately, and then combined with the raster through some scripting- but that seems messy.

If this is not possible to do with geoserver alone - perhaps there is an extension that would help? If not, then could someone recommend another type of web service that could do such a thing?


texts: Adding a popup to features in a leaflet.vectorgrid layer

texts:
I am loading a tile set from MapTiler to my leaflet map. I have no trouble styling each feature in the tile set. I can set each feature to the appropriate color.

I would like to also add a popup to each of the features. Is this possible? How? Where is the documentation that explains this?

const FOO_TILES_PBF_URL =  "https://api.maptiler.com/tiles/<tileset_identifier>/{z}/{x}/{y}.pbf?key=<my_maptiler_key>";

  function getColorFromRedBlueRange(d) {
    return d == 8
      ? {
          fillColor: "yellow",
          fillOpacity: 0.7,
          color: "yellow",
          fill: true,
          weight: 1.5
        }
      : d == 7
      ? {
          fillColor: "red",
          fillOpacity: 0.7,
          color: "red",
          fill: true,
          weight: 1.5
        }
      : {
          fillColor: "black",
          fillOpacity: 0.7,
          color: "black",
          fill: true,
          weight: 1.5
        };

var FooTileOptions= {
    interactive: true,
    vectorTileLayerStyles: {
      LayerName: properties => {
        return getColorFromRedBlueRange(
          properties.foo_variable]
        );
      }
    }
  };

  var fooLayer= L.vectorGrid.protobuf(
    FOO_TILES_PBF_URL,
    FooTileOptions
  );


texts: Creating a dissolved buffer layer in PostgreSQL/PostGIS

texts:
I have got a MultiPolygons layer to which I want to create 10 meter buffer around and then dissolve. All data are having SRID 4326.

I am new to PostgreSQL/PostGIS. I am drifting around the ST_Buffer and ST_Union commands, but not reaching to any conclusion.

Here is my code:

CREATE TABLE public.Buffer10mtr AS
  SELECT ST_Union(ST_Buffer(Zone2.shape as geometry,10)as geometry)
  FROM public.Zone2;

What can I change?


texts: What are the requirements for using Google's GEvent?

texts:
For a current project, we need to know when a Google map is done loading all of its tiles. I found the following question, which was for the same basic problem: register event "loadend" on layer

So now, we're using the code from that answer:

GEvent.addListener(mGoogleHybridLayer.mapObject, "tilesloaded", function() {
 console.log("Tiles loaded");
});

However, this gives a reference error for GEvent. Are there some extra requirements for using GEvent? Do we need to import it from somewhere?

I can't post code we're using, because there isn't any yet - even in the browser's debug console, GEvent is undefined.

Using Openlayers 2.11.

EDIT: maybe I should restate the question into one that can actually be answered: what are all the requirements for being able to use GEvent?


texts: Calculating continuous raster representing land cover proportion from categorical raster using QGIS?

texts:
I have a raster (5 x 5 m resolution) representing land cover. Each cell contains a integer value (from 1 to 10) representing a total of ten land cover classes (e.g. 1-forest; 2-agricultural areas; 3-urban areas;(...)10-scrublands).

Based on this categorical land cover raster, I would like to produce ten rasters with a resolution of 10 x 10 m representing the proportion of each land cover class per cell. In another words, I would like to produce a raster representing the proportion of forest, another raster representing the proportion of agricultural areas, etc.

This figure explains exactly what I want to do (for illustrative purposes, I only accounted for 3 classes in the figure).

I have been looking for similar questions in the internet, but usually people had slightly different questions (e.g. calculate land cover proportion within a polygon).

Does anyone know how to do this in QGIS?


texts: Get grid area from geographic coordinates

texts:
I have a list of geo coordinates. What I want is a list of grid areas and a number of geographic points in that area. By doing this I can get the information in which areas in the world the most points are set.

Example:

[lat, lng], ... [lat, lng] -> [gridId, total number of geo coords in that grid]

I want the size of the grid to be configurable, e. g. Zoom Level 8, or by saying Grid with/height should be e.g. 10km. Since I can have a huge amount of positions (possibly hundredthousands of points) the calculation should be efficient. Ideally the list of [gridId, totalNr] should be sorted by totalNr.

What I am currently doing is to use DBSCAN clustering to cluster geo points into groups and calculate the centroid of that group. I think that the grid algorithm could be faster, if the calculation of grids is fast enough.

Im doing it in Javascript. I wanted to use elastic search, but could not find appropriate functions. I also could not find other projects for such an use case.


texts: Aligning existing polygon borders using QGIS?

texts:
I'm looking for a method to align the borders of two overlapping and already existing polygon layers. I'm aware of the node tool and its tolerance snapping option, but the area I'm working with is too large to do it manually.

I have two polygon layers. Layer 1 contains about 11,000 municipality features. Layer 2 contains 14 climatic zones (colored in picture 1).

As seen in picture 2 the outer border of the two layers do not fit. First of all I would like to manipulate the outer border of the climate zone layer (colored) so that it aligns with the outer border of the municipality features. Further I do not want the climate zone features to split the municipality features. One municipality shall be covered by only one climate zone.

At the very end I plan to intersect the two layers so that I have a list: Municipality ID to climate zone ID.

So far I only found a tool for ArcGIS, called Align to Shape.

Is something like this available for QGIS too?


texts: Sending carto a SQL API request for column names returns error

texts:
I got a tip from a previous question for getting a list of column names in a Carto table.

I enter this into my browser: https://MY_CARTO_USERNAME.carto.com/api/v2/sql?q=SELECT column_name FROM information_schema.columns WHERE table_schema = 'MY_CARTO_USERNAME' AND table_name = 'TABLE_NAME'&api_key=MY_API_KEY. Of course, I replace MY_CARTO_USERNAME and MY_API_KEY with my Carto username and API key.

This is the error I get: {"error":["You must indicate a sql query"]}

What must I change in my query to get the column names of my Carto table?


texts: Mouse right-click is not working anymore when using ArcGIS10

texts:
The mouse right-click option has stopped working with ArcGIS10. It happened suddenly - the mouse is ok (not broken) as it works with other programs. Do you know what happened and how to fix it?

The right-click option is useful for editing as you can just right-click select from the drop down list, absolute coordinates, direction, delete sketch etc. It is also not working, when I'm in layout view as I can't change grid properties.


texts: Convert a vector file to an RGB raster using gdal_rasterize without filling polygons

texts:
I am trying to convert a vector file that includes polygons to an RGB raster using gdal_rasterize without filling the polygons.  Is there a way to do this?

For example, I use the following command to create an RGB raster with opaque black values showing the vector information from a vector file:

gdal_rasterize.exe -burn 0 -burn 0 -burn 0 -burn 255 -ot Byte -ts 927 583 -te -87.7595200000000 41.8042200000000 -87.5067400000000 41.9631900000000 -of GTiff -l ACHARE -l BCNLAT -l BCNSPP -l BERTHS -l BRIDGE -l BUISGL -l BUAARE -l BOYLAT -l BOYSPP -l CBLARE -l CBLOHD -l CBLSUB -l CTNARE US5IL11M.000 out.tif

When I run the above command, I end up with an RGBA raster with a single black box covering most of the image because the vector file has a polygon bounding all the other layers.  Is there a way to rasterize only the boundary of this polygon instead of all the area inside this polygon?  Is there another tool (preferably with GDAL) I could use to do this?

I might be able to get all the point and line layers in my vector file and only display these layers.  This is one possible solution for me if there's not an easy way to do this using gdal_rasterize.  I might also be able to convert the polygons to lines and rasterize them as lines.  That's another possibility.  Both of these solutions would be more work than I was hoping to do to solve this problem.


texts: How to save points without losing z dimension

texts:
I try to save Points with OGR. It works well but I lose the z dimension while saving/loading.

from osgeo import ogr
from osgeo import osr

shapePath="...xxx.shp"
X='1.0'
Y='2.0'
Z='3.0'

driver = ogr.GetDriverByName('ESRI Shapefile')

srs = osr.SpatialReference()
srs.ImportFromWkt(mySRSinwkt)

dataSource = driver.Open(shapePath, 1)
layer = dataSource.CreateLayer("points", srs, geom_type=ogr.wkbPoint)
layer.CreateField(ogr.FieldDefn("Point_ID", ogr.OFTInteger))

# Set feature fields
feature = ogr.Feature(olayer.GetLayerDefn())
feature.SetField("Point_ID", 1)

# Set geometry
geom_wkt = 'POINT (' + X + ' ' + Y + ' ' + Z + ')'
point = ogr.CreateGeometryFromWkt(geom_wkt)
feature.SetGeometry(point)

#Point has z dimension
print str(feature.GetGeometryRef())

# Create feature
layer.CreateFeature(feature)

# remove feature and datasource from memory
feature.Destroy()
dataSource.Destroy()

# Load
dataSource = driver.Open(shapePath, 0)
layer = dataSource.GetLayer()

# Point has no z dimension anymore
for feature in layer:
    print str(feature.GetGeometryRef())

dataSource = None


texts: Preparing geometry for psycopg2 using ArcPy without SHAPE@WKT?

texts:
Trying to print out a polygon shapefile in WKT format in a standalone script. I am using ArcGIS 10.1 and use the da modules all the time. I have seen this post Converting geometry to WKT using ArcPy? and have copied and pasted that code to use on my polygon shapefile and it does not work. The error I get is

RuntimeError: unknown geometry token 'SHAPE@WKT'

import arcpy
arcpy.env.workspace = "C:\\Users\\Ralph\\Desktop\\sql"
shp = "C:\\Users\\Ralph\\Desktop\\sql\\states2.shp"
for row in arcpy.da.SearchCursor(shp, ["SHAPE@WKT"]):
    print row[0]

I know it is mentioned in the above post that this token was made available for service pack 1. I am starting to think I do not have service pack 1 nor do I know how to check if I have it. I plan on reading this polygon in a easy format so I can import it using psycopg2 to postgreSQL. Is there another way to read polygon geometries for this task only by using Arcpy?


texts: How to store / what algorithm to use to find closest points from a given set of points to a polyline layer

texts:
I have a polyline layer and a points layer with features that are nearby the polylines. For the latter I want to find the closest points on the polylines (as existing vertices or new-created vertices along the polylines).

My initial approach is to loop through all polylines and check if their bounding boxes with a bit of buffer are intersecting the given points, so 2 for-loops then for each polyline candidate to loop through vertices and apply the distance point-line algorithm to find the closest point. But can become very slow especially for large datasets.

I've quickly checked Quadtrees but I see they are optimized for storing and fast lookup of point features. Is there anything similar for polylines or rectangles (bounding boxes)?

Later edit: I'm more interested in the algorithm / ways to store data because it's part of the software that I am implementing at this moment. So any existing GIS software / plugins examples are ok but I want to see what's "under the hood".


texts: How to extract a river profile using specific distance intervals

texts: I'm using QGIS 2.16 and I'd like to extract a river profile from a DEM applying user-defined intervals (e.g. every 100 m of horizontal distance). I get profiles using both the "Profile tool" and the "Profile from lines" SAGA process. However, I cannot control how those profiles are generated.


texts: Incorrect Polygon Area Calculations in ArcGIS Desktop?

texts:
I have converted a raster to a polygon, dissolved the polygons into singlepart, and tried to calculate areas using Calculate Geometry. Many of the areas are 'zero' or are areas that are obviously not correct.

How can I fix this?

When I calculate the areas for the undissolved polygons they appear to be correct. They are just subdivided into grid cells with the same spatial resolution as the Landsat image from which they are derived (30m x 30m). So, one grid cell has an area of 900 sq meters. The problem is that I need the grid cells dissolved so I can determine the exact number of polygons and their respective areas, hence the use of the dissolve tool.

The defined projection is WGS UTM Zone 9N (the same as the Landsat image) and it has never been changed.


texts: Identify spectral bands in a multispectral raster dataset

texts:
I'm trying to make an application to use on precision agriculture using gdal, and want to give it some flexibility to support imagery from different providers.

I've been looking around if there's a way to identify the bands in a multispectral image despite the provider.
I used gdalinfo with some raster datasets from Landsat, Spot6 and Sentinel-2, but the only hint it gives I could use is Color interpretation, which is not very useful.

I get gray for the individual bands of the Sentinel-2 image, and Red, Green, Blue, and Unidentified for landsat and spot6 stacked datasets I prepared on ArcGIS.

Is there a way to identify each band from the metadata, the header of the tiff image or something, so my program can identify each band without me telling it each time I use a dataset from a different provider?


texts: How important is a knowledge of computer science in GIS analysis?

texts:
Don’t think this question has been asked before; I’ve just come across a question on how much math one has to know for a career in GIS. As a geography student I’m wondering how much computer science I’d need to know if I’m planning on being a GIS analyst, particularly as nearly any answer I read on here quickly delves into how to solve something through programming and it’s beyond me. Furthermore, I’ve noticed both on here and offline that most get into the GIS line of work through computer science rather than geography. I’m nearing the end of my undergrad and will be doing an MSc in GIS and remote sensing but I have also been told by a career advisor to look into a conversion course in computer science before I do the Masters to give me an edge.

However, I am apprehensive of a few things, the main reasons for my qualms being:

Aptitude.
I’m wondering if I, as just a French/geography student, am cut out for computer science. More specifically, how much math is involved? While I am not averse to it (and have been using light statistics and quantitative methods increasingly in geography for my final year) it’s probably not my strongest suite. What other skills are crucial? Does computer science require more abstract of concrete thinking, as I think I might slightly lean towards the former?

The actual point?
Though it’s a more applied Masters, I’ll be doing a little programming for the course. Would doing an entire year of computer science be an addition to me, not only in preparation for the Masters but also long-term, keeping in mind that it’s not only more fees but also will take a lot of effort? It’s known to be one of the most intense post-grads in the university. It is, after all, teaching the bulk of a three-year Bachelor’s in one year. I’m wondering if it will pay off. If it would be beneficial, I'm willing to spend an extra year in education.

To give an idea of what the course will entail, here’s the list of modules:
•   Introduction to computer systems (both semester 1 and 2)
•   Algorithms and data structures (both semester 1 and 2)
•   Databases
•   Structured programming
•   Discrete structures
•   Web information processing
•   Software testing
•   Multimedia technology
•   Software engineering and software processing
•   Operating systems


texts:
Python was working and then suddenly

I have attempted to find /Users//Library/Preferences/org.qgis.QGIS2.plist as advised in Failed to open Python console error message but this is not in my preferences folder.
I am running Las Palmas 2.18.14 on Mac OS High Sierra 10.13.2.

Is there a way to troubleshoot the problem or does anyone have any ideas on how to fix this?


texts: Placing image inside of Label using ArcGIS Desktop?

texts:
In ArcView, can I add an image to a label?

Assume I have the photograph path as an attribute.


texts:
I have some code that looks something like this:

def get_objs(all_objs, line, my_distance):
    """
    Gets all objs within my_distance of line.

    all_objs is a numpy array of shape (x, 4) denoting bounding-box coordinates.
    line is a LineString object with an arbitrary number of points.
    """
    line_objs = []

    for o in all_objs:
        x_avg, y_avg = np.mean([o[0], o[1]]), np.mean([o[2], o[3]])

        p = Point(x_avg, y_avg)  # Point is a Shapely type

        if p.distance(line) <= my_distance:
            line_objs.append(o)

    return line_objs

My question is whether I can broadcast in this case, because this code is really slow on my data and I suspect a vectorised broadcast would help a ton. I have checked the numpy docs, but they don't make it clear that you can broadcast your inputs to a new type and then do something with that type.

Likewise the Shapely docs don't make it clear that you can apply the .difference() function in element-wise fashion so that you can compare to some geometric object. Perhaps I am being a bit optimistic. What I have in mind is something like:

x_avgs = (all_objs[:, 0] + all_objs[:, 1]) / 2
y_avgs = (all_objs[:, 2] + all_objs[:, 3]) / 2
xy_avgs = np.array((x_avgs, y_avgs)).transpose()

# element-wise initialisation and comparison, followed by
# masking for true-y elements.
line_objs = (Point(xy_avgs).distance(line) < my_distance).nonzero()

return line_objs

The above does not work, since Point()s don't accept arrays :(

After some (very statistically shaky) benchmarking, a quick consultation with timeit:

>>> setup = "import numpy as np; a = [(x, y, x, y) for x in range(100) for y in range(100)];"
>>> print(timeit.timeit("[(np.mean([o[0], o[1]]), np.mean([o[2], o[3]])) for o in a]", setup=setup, number=100))
30.91517815797124
>>> print(timeit.timeit("[((o[0] + o[1]) / 2, (o[2] + o[3]) / 2) for o in a]", setup=setup, number=100))
0.3588130779680796

... That's insane! That's nearly 2 orders of magnitude on a sample much smaller than mine. Perhaps vectorisation isn't necessary after all.


texts: Transfer multiple point to polyline based on spatial location

texts:
I have a number of polylines (bus routes) and points (bus stops). The polylines span from stop to stop. Each stop has an ID, a single number.

I am trying to get the point ID that touches each polyline into the polyline so each polyine will have stop 1 and stop 2 which it spans between.

Any thoughts on how I can do this?

I have tried a spatial join which will give me two polylines, one for each stop it intersects.


texts: Location issues in solving OD Cost Matrix ArcGIS 10.6.1

texts:
I created a network dataset only for cars from my overall street network.
I am getting a partial solve for The two warning messages I keep getting are Warning: Location "0935010000038273" in "Origins" does not have a corresponding network element. The network dataset needs to be rebuilt.

and also that there is no destination for location "x." I have 61814 point origin locations (all of which display as located, no errors) and 1 point destination. The OD cost matrix solves for 53213 locations out of the 61814. There are no barriers applied, and the only restriction I have is the Oneway restriction automatically generated by the Network Dataset wizard.

The parameters I used for adding the origin and destination locations was to use the network location information for each location, and I set all point location data to snap to the network with a 5m offset.

Is this an issue with the network dataset restriction or the network dataset itself? Or is it an issue I did not have the error message about corresponding network elements for using the same network (but selected for cyclists, and no one-way restriction). I previously ran an OD cost matrix for cycling that solved for majority of locations, but had also the "no destination for location x" warning.


texts: How do you remove the title in the banner from the ESRI geoportal extension?

texts:
I would like to remove the catalog.site.title text from the banner area in the geoportal page. However, I would still like to keep the title of the page in the browser bar/window.

Does anyone know in which file the catalog.site.title is located for just the banner text? If I remove the text from the catalog.site.title from the site resources in the gpt.properties file the banner title and browser title are replaced with
???catalog.site.title.


texts: Vector Grid shows big difference between cartesian and ellipsoid measurement

texts:
I am using QGIS 3.4.1 on Windows 10. Project CRS is EPSG:3857 (WGS 84/Pseudo-Mercator).

I have defined a vector grid. Type is Rectangle (polygon). Horizontal and vertical distance is 50 m. Overlay is 0 m. GRID CRS is EPSG:3857. This goes well. But after selecting one of the polygons, the Cartesian properties are quite different from the Ellipsoid properties:
- Cartesian: Outline: 200 m, Area 2500 m2
- Ellipsoid: Outline: 122,536, Area 938,435 m2
- Using the measurement tool shows sides of the rectangle of 30,567 m.

I expected a small difference between Cartesian and Ellipsoid outcomes, but this is impossible.

Does anybody have a clue of what is going wrong?


texts: Querying OpenStreetMap tags by polygon using QGIS?

texts:
If I understand OpenStreetMap (OSM) correctly, all objects essentially have associated tags based on type of object (street, building etc) ...

Can I use QGIS to obtain the statistics of different commonly used tags in that area so that I can have an idea of the feature distribution?


texts: Coordinate transformation differences in ArcGIS, QGIS and Google Earth Pro

texts:
I have a shapefile with points in Adindan UTM 37N coordinates. Now there seem to be roughly two ways of converting those coordinates to WGS84.

When loaded in Google Earth Pro the points are about 150m shifted to the southwest compared to in QGIS.

In ArcGIS both of these behaviours are possible. When using on-the-fly projection with 'Adindan_To_WGS_1984_1' selected in the Data Frame Coordinate System Transformation the result is the same as in QGIS, (both using towgs84 geocentric translation -166,-15,204).

When <None> is selected however the result is the same as in Google Earth. (This is only possible in on-the-fly projection, the Project tool always selects the closest match by default, just as with KML export.)

This Google Earth/ArcGIS <None> behaviour is a black box for me, but in this case of converting Garmin GPX data in Adindan it seems to be the correct behaviour, so my question is: What happens in this case, and is it possible to do the same transformation with QGIS/proj4?
(I tried using +towgs84=0,0,0,0,0,0,0 but that causes an even bigger shift in a different direction)

For example: The Adindan coordinate [391000, 1221000] is converted to [38.002975, 11.045626] by QGIS and to [38.002148, 11.044758] by Google Earth.

Another way to use the coordinates is to interpret them as being WGS84 UTM 37N, then they transform to [38.002132, 11.043761]


texts: Select GeoTIFF from a database via GeoServer and OpenLayers3

texts:
I'm almost a beginner at PostGIS and GeoServer so I hope you could help me with my question.

In my PostGIS database I have a lot of GeoTIFFs that I would like to display in a WebGIS via GeoServer and OpenLayers3.

Every hour a new GeoTIFF with temperature information is coming.

INCA-T2M_20160614T1300.tif

INCA-T2M_20160614T1400.tif

INCA-T2M_20160614T1500.tif etc.

I upload them via raster2pgsql into the database.
There they have the columns rid, rast and filename

Now I would like to build a query in the WebGIS where the user could type a date and time and the suitable GeoTIFF will be displayed.
Therefore I have to build the correct WMS request with OpenLayers3 and send them to the GeoServer.
But how does the correct WMS-request looks like and how the GeoServer will ask my database for the correct table?

After reading the GeoServer documentation maybe the creation of an sql view with parameter could be the solution?

The examples here are dealing with the selection of a specific range of data from a fix table. But I don't know how do I have to build a query to ask for a specific table, not for a range of data in a fix table.

But maybe I'm on a wrong path... Do you have any idea?


texts: Changing the icon picture on a custom toolbar button

texts:
I've made a few customised plugins using Plugin Builder and QT. They're all up and running, but I'm having some trouble changing the icon picture on the buttons that I've added to the toolbars to run my plugins.

Is there a way to insert my own picture in the button, or even have text instead? The only pic that seems to work is the default icon.png created by the Plugin Builder. I've tried changing the path in the python script (see below) but no avail, it only works with icon.png. Does a new picture have to be a certain size or type?

def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    icon_path = ':/plugins/button/icon.png'
    self.add_action(
        icon_path,
        text=self.tr(u'button'),
        callback=self.run,
        parent=self.iface.mainWindow())


texts: ArcGIS Desktop 10.2.1 and one-to-many relationship

texts:
I've been struggling to link a feature of 450 municipalities and approximately 60,000 fires that have occurred over the last 15 years. I've attempted to make a relationship class. I can see the individual events when I query a municipality but I cannot symbolize based on, perhaps, graduated symbols to represent the number of fire events normalized by size. I've attached a pic of the tables and the Identify.

Not sure if I could summarize the number of instances of fires and append them as an attribute to the municipal boundaries instead..? I'm in an FGDB, no DB backend.


texts: Automatic split a mosaic of satellite images into pieces

texts:
I have a mosaic of satellite images, but for the analysis, I need to define the boundaries of each part. Does anyone know how to do it in an automatic way?
I need to get the red boundaries like on the second image but in an automatic way(please do not offer a manual digitization option).

I'm looking for some tool or algorithm that can help to make those boundaries automatically (there's a lot of similar). Maybe some raster processing techniques or should I dig in machine learning?


texts: Not able to load mbtiles in ArcGIS runtime SDK

texts: I'm not able to load mbtiles in any layer using arcgis runtime sdk version 100.1 for iOS. Is it possible to do that in new version of run time sdk. In old runtime SDK, I have used AGSTiledServiceLayer to load mbtiles in version 10.2 and now can load only tpk files using AGSArcGISTiledLayer.


texts: Displaying raster data export from temp folder using ArcGIS Desktop?

texts:
I'm coming up against a fairly unusual raster display issue I haven't encountered before. As a caveat, I'm not the most experienced at working with display of raster data, so I assume (and hope) it's a simple solution.
I took a DEM and created a hillshade from it using the image analysis window, so far so good, made some adjustments to transparency & contrast. This was stored in my temp folder. The result:

When I exported to the folder I am keeping all my data, however, (right click temp layer in TOC etc., used default settings in export window) and added the resultant TIFF to the map, it looked like this:
Note that the original layer is still turned on- the only output was the black bands. I've calculated statistics for the layer hoping for a quick fix but no luck. The histogram seems to indicate there's some data there other than just black pixels in seemingly arbitrary strips.

But there's definitely a clear trend in the data bunching up at the 180 mark.
Is there something incredibly obvious I'm missing? I'm not sure if this is a display issue or something unexpected happening to the data during the export process.
Using: ArcGIS 10.3 on Windows


texts:
I've just tried some simple code from example:

$(document).ready(function() {

   var map = OpenLayers.Map('map-content', {});
   var wms_serever_map = OpenLayers.Layer.WMS('Base l', 'http://vmap0.tiles.osgeo.org/wms/vmap0', {layer:basic},
       {});
   map.addLayer(wms_serever_map);
   if(!map.getCenter()){
      map.zoomToMaxExtent();
   }
});

but chrome and firefox don't show me the map and throws error: uncaught type error: Cannot call method 'apply' of undefined Openlayers.js:127 Its rather interesting, because some time ago this code worked


texts: Serving geospatial files from a file system using GeoWebCache

texts:
Is there a primer somewhere for how to go from having various geospatial files (such as shp, dbf, gdp, geotoff, geojpg, geopng) to being able to serve them to a map via service calls to GeoWebCache directly from a file system without having to put something like GeoServer in the loop?

Assuming this is even possible, and I believe it is, but it's not easy to figure out how to do this from the documentation.

Can someone confirm whether or not this is possible, and if so list the steps that would be needed?


texts: Using gdalwarp to merge bands into a mosaic created by gdal_merge?

texts:
I have created a mosaic with gdal_merge.py using the following command:

gdal_merge.py -createonly -separate -co PHOTOMETRIC=RGB -o merged.tif LE07_L1TP_046030_20180810_20180810_01_RT_b3.tif LE07_L1TP_046030_20180810_20180810_01_RT_b2.tif LE07_L1TP_046030_20180810_20180810_01_RT_b1.tif

It seems to create the mosaic no problem.  Now, when I run the gdalwarp command to create a viewable image (.tif), I get a red image, even though I am merging the color bands.  I am using this command:

gdalwarp --config GDAL_CACHEMAX 512 -wm 4096 -nomd LE07_L1TP_046030_20180810_20180810_01_RT_b3.tif LE07_L1TP_046030_20180810_20180810_01_RT_b2.tif  LE07_L1TP_046030_20180810_20180810_01_RT_b1.tif  merged.tif

The output of gdalinfo on merged.tif only gives me a red color band, where the other color band values are 0.  Does anyone have any ideas?  Also, I know I can use gdal_merge to do all of this at once, but I have a memory constraint that prevents me from using gdal_merge to do all of it at once.

gdalinfo output:

    Size is 9930, 9315
Coordinate System is:
PROJCS["Albers",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Albers_Conic_Equal_Area"],
    PARAMETER["standard_parallel_1",29.5],
    PARAMETER["standard_parallel_2",45.5],
    PARAMETER["latitude_of_center",23],
    PARAMETER["longitude_of_center",-96],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-2336089.094018519900000,2701776.184199890100000)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-2336089.094, 2701776.184) (125d42'39.19"W, 43d59'36.72"N)
Lower Left  (-2336089.094, 2422326.184) (124d40'18.81"W, 41d36'34.34"N)
Upper Right (-2038189.094, 2701776.184) (122d 7'21.46"W, 44d46' 3.96"N)
Lower Right (-2038189.094, 2422326.184) (121d11'45.22"W, 42d21'15.67"N)
Center      (-2187139.094, 2562051.184) (123d25'35.42"W, 43d11'34.14"N)
Band 1 Block=9930x1 Type=Byte, ColorInterp=Red
    Computed Min/Max=0.000,255.000
  Minimum=0.000, Maximum=255.000, Mean=27.268, StdDev=43.782
  Metadata:
    STATISTICS_MAXIMUM=255
    STATISTICS_MEAN=27.267951614063
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=43.781663272236
Band 2 Block=9930x1 Type=Byte, ColorInterp=Green
    Computed Min/Max=0.000,0.000
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0
Band 3 Block=9930x1 Type=Byte, ColorInterp=Blue
    Computed Min/Max=0.000,0.000
  Minimum=0.000, Maximum=0.000, Mean=0.000, StdDev=0.000
  Metadata:
    STATISTICS_MAXIMUM=0
    STATISTICS_MEAN=0
    STATISTICS_MINIMUM=0
    STATISTICS_STDDEV=0

Is it overwriting the red band each time?  If so, is there a way I can select which color band each b*.tif goes to?


texts: Why do I get "Insufficient memory for image" although there should be enough RAM?

texts:
QGIS reports "Insufficient memory for image 1403x716", when I want to load a raster. The problem begins after 1 hour of work.

Had somebody had such problem?

System  win7 64 bit, 16GB RAM, Intel Core i5.


texts: GDAL raster shifting with python error: 'NoneType' object has no attribute 'GetGeoTransform'

texts:
I'm trying to shift a georeferenced raster by 0.5m N and 0.5m W in QGIS using the python console:

from osgeo import gdal

# Open in read/write mode
rast_src = gdal.Open('test_raster.tif', 1)

# Get affine transform coefficients
gt = rast_src.GetGeoTransform()

# Convert tuple to list, so we can modify it
gtl = list(gt)
gtl[0] -= 0.5  # Move west 0.5 m
gtl[3] += 0.5 # Move north 0.5 m

# Save the geotransform to the raster
rast_src.SetGeoTransform(tuple(gtl))
rast_src = None # equivalent to save/close

However line #7 (gt = rast_src.GetGeoTransform()) throws an AttributeError: 'NoneType' object has no attribute 'GetGeoTransform'. The installed GDAL version is 1.11.0. The raster seems fine. So what's wrong?


texts: Specifying input coordinate reference system for WFS?

texts:
GeoServer allows to specify output CRS by parameter srsName in WFS GetFeature.

Is there a way to specify the input CRS?

I am displaying layers of different CRSs in my map - all reprojected to EPSG:900913. Therefore the WFS requests generated are also in EPSG:900913. What I want is to tell GeoServer that my inputs are in EPSG:900913 so that GeoServer can reproject my input and return the result.

I am using WFS version 2.0.0.


texts: python GDAL raster calculation error message

texts:
hello I try to follow this post and Friend commend :

import os
gdal_calc =  'python C:/PROGRA~2/GDAL/gdal_calc.py -A d:/mnt_5x5.tif --outfile=d:/test.tif --calc="A+1"'
os.system(gdal_calc)

but I take error :

..<string>:1: RuntimeWarning: overflow encountered in multiply
    C:/OSGeo4W64/bin/gdal_calc.py:296: RuntimeWarning: invalid value encountered in multiply

here complete message :

0 ..<string>:1: RuntimeWarning: overflow encountered in multiply
C:/OSGeo4W64/bin/gdal_calc.py:296: RuntimeWarning: invalid value encountered in multiply
  myResult = ((1*(myNDVs==0))*myResult) + (myOutNDV*myNDVs)
 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 .. 90 .. 100 - Done

any idea why I take this error ?I take new image for output but I am not sure if this image is correct after this error .

here my complete code :

gdal_calc =  'python C:/OSGeo4W64/bin/gdal_calc.py -A {0} -B {1} -C {2} --outfile=test.tif --calc="(5*A)+(5*B)+(2*C)"'.format(a,b,c)
os.system(gdal_calc)


texts: Convert vector data to raster data in GeoTools

texts:
In my project I am able to open and save shapefile. I want to convert this vector data to image format as a TIFF image.

Please help me the procedure using GeoTools.


texts: Color Correction from two nearby Landsat8 Rasters in QGIS?

texts:
I have two Landsat8 images and they are side by side, but they are from different Dates. They look similar, but the color of the same forest is a bit different.

I tried a color correction in QGIS with the tool "r.blend - Blends color components of two raster maps by a given ratio" from GRASS GIS.
The result is not what I expected and wanted and the tool even combine the rectangular frame from the Raster.

Is there any other tool in QGis for a color correction from two Rasters?

I am a beginner in QGIS and I am not a native English speaker.


texts: Converting Vertices of Polylines and Polygons to Points and Maintain Attributes

texts:
I'm pretty sure someone has done this before, but I'm having no luck.  I'm using ArcGIS Standard version 10.2.1 and I want to convert all the vertices of a polyline or polygon to points.  Essentially, I want to emulate the Feature Vertices to Points tool (which I can't use due to my ArcGIS license).  I'm writing a python tool, so any answers should reflect a python/arcpy solution.  Therefore, nothing like ET Geowizards, etc.

I want to maintain the fields of the original data and in addition, add X and Y fields as well as a point identifier.  So, if the input fields look like this:

LABEL |  TYPE  |  CODE
abc   |  type1 |  3
def   |  type2 |  4

Then the resulting output points might look like this:

LABEL |  TYPE  |  CODE  |  LONG  |  LAT  | ID
abc   |  type1 |  3     |  153   | -27.5 | 1
abc   |  type1 |  3     |  150   | -27.5 | 2
def   |  type2 |  4     |  151   | -30   | 1
def   |  type2 |  4     |  151   | -29.5 | 2
def   |  type2 |  4     |  151   | -30.7 | 3
def   |  type2 |  4     |  151   | -30   | 4

I've looked at this solution, which is fantastic and quick, but it doesn't maintain the fields of the original data.  I've also looked at modifying this blog post which seems to be the most promising, but again, I can't figure out how to keep the original fields.

My code looks something like this:

fields = arcpy.ListFields(input_fc)
arcpy.CreateFeatureclass_management(destination_workspace, point_file_name, "POINT")
for field in fields:
    fieldnames.append(field.name)
    if field.type in ('TEXT', 'FLOAT', 'DOUBLE', 'SHORT', 'LONG', 'DATE', 'BLOB', 'RASTER', 'GUID'):
        arcpy.AddField_management(point_file_name, field.name, field.type, field.precision, field.scale, field.length, field.aliasName, field.isNullable)

fieldnames.append("SHAPE@JSON")
cursor = arcpy.da.InsertCursor(output_fc, fieldnames)

with arcpy.da.SearchCursor(input_fc, fieldnames) as c:
        for row in c:
            geom = json.loads(row[-1])
            for path in geom['paths']:
                for pt in path:
                    cursor.insertRow(row[0:-1],pt)

It's failing with the insert cursor.  The error message I'm getting is:

  Traceback (most recent call last):
  File "H:\Python\CreatePoints_v0_2.py", line 60, in
  cursor.insertRow(row[0:-1],pt) RuntimeError: Cannot find field 'LABEL'

The input feature class has a field called "LABEL".  However, it's not adding that label to the created feature class.

Note, I want to convert vertices, not centroids.

I'm probably going about this the wrong way, so if anyone has ideas on how to either modify my code, or point me in another direction, I'd be appreciative.


texts: Interpolating species data from sample points in a small area

texts:
I have sample data from a small site (10 acres) which shows the number of occurrences of a particular species of snail over two sampling periods (fall and spring sampling). I want to try and interpolate the potential concentration across the entire site based on the sample data.

Any suggestions on how to best do this?

IDW?

Kriging?

Density?


texts: Exporting time series data of multipoints in Google Earth Engine

texts:
I have created three points of land surface temperature. I want to export the time series result in three different variable names.
// Import country boundaries feature collection.

var dataset = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017');

// here is the three points of land surface temperature

var features = [
ee.Feature(ee.Geometry.Point(14.5, -15)),
ee.Feature(ee.Geometry.Point(14.8, -15.1)),
ee.Feature(ee.Geometry.Point(14.9, -15.2))
];

var points = ee.FeatureCollection(features);
Map.addLayer(points);

// Import LST image collection.

var modis = ee.ImageCollection('MODIS/MOD11A2');
var start = ee.Date('2001-01-01');
var dateRange = ee.DateRange(start, start.advance(15, 'year'));
var mod11a2 = modis.filterDate(dateRange);

// Select only the 1km day LST data band.

var modLSTday = mod11a2.select('LST_Day_1km');
var modLSTc = modLSTday.map(function(img) {
return img
.multiply(0.02)
.subtract(273.15)
.copyProperties(img, ['system:time_start']);
});

// Time series at the identified points

// "Please help here in this point"

The export result looks like the attached picture.


texts: Combining different svg fill symbols in QGIS without overlapping

texts:
Is it possible to get a polygon filled with a couple of different SVG fills without overlapping?
I want them to go side by side forming a mosaic, instead of overlapping each other, but can't find how (of if) I can do that. I guess there should be something like an offset option, but there's no such thing on SVG fill tab, or I don't know how to find it.


texts: Creating a table that stores Distances and other column

texts:
I am looking to run a query to calculate the distance from a polygon to other features (will use this query on point, polygon, and line features) within 2 miles and create a table that will display multiple columns:

Polygonid
nearfeature_dist
nearfeature_name

Currently I have:

CREATE TABLE near_fiber AS
SELECT parcels.objectid,

Not sure what postgis function I should use to get my intended results. Im thinking a mix between ST_Distance and ST_DWithin.

I will have to use this query multiple time for different features.


texts: Creating grid of latitudes and longitudes to be displayed in polygon boundary with ArcMap?

texts:
I have with me a polygons of a country, state and city in geographic coordinate system. A grid of latitudes and longitudes is to be created to be displayed within the polygons. The interval between the graticule grid lines will vary depending on the extent (country/state/city) for which the grid is to be created. The grid needs to be a shapefile or feature class since a conversion to kml/kmz is later required.

Are there any tools in ArcMap that will allow to achieve this?


texts: Importing a georeferenced DWG/ DXF into QGIS

texts:
I work in an ecological consultancy and we would save money on mapping tiles if we could import DWG files into QGIS, they must be georeferenced however as we need to draw on various habitats accurately and be able to measure distances/ areas.

I have obviously managed to convert DWGs into DWF using Thienga and opening it into QGIS, but it loses the georefrence data.

I have also tried Project > import DWG/ DXF - but cannot find the 'geopackage' (whatever that is)?


texts: Using IAddDataDialog2 to obtain layer from .gdb with ArcObjects and C#

texts:
I am using ArcMap 10.4.1.

I cannot figure out how to properly initialize/create IAddDataDialog2 in C# to access a .gdb file in which I'm suposed to find a layer that I want to add in my TOC. I tried reading out the API from:

IAddDataDialog2
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#/IAddDataDialog2_Interface/000v00000025000000/

IAddDataDialog
http://help.arcgis.com/en/sdk/10.0/java_ao_adf/api/arcobjects/com/esri/arcgis/arcmapui/IAddDataDialog.html

AddDataDialog
http://help.arcgis.com/en/sdk/10.0/java_ao_adf/api/arcobjects/com/esri/arcgis/arcmapui/AddDataDialog.html

But I still cannt figure it out.

EDIT: PolyGeo's request

 /// <summary>
 /// Do stuff.
 /// </summary>
 /// <param name="sender">The sender.</param>
 /// <param name="eventArgs">The <see cref="EventArgs"/> instance containing the event data.</param>
 private void BoutonAjouterCoucheOnClique(object sender, EventArgs eventArgs)
 {
      ServiceAddDataDialog service = new ServiceAddDataDialog();
      service.OuvertureDuDialogAddData();
 }

/// <summary>
/// Service de support pour la fenêtre d'ajout de données.
/// </summary>
public class ServiceAddDataDialog
{
    /// <summary>
    /// Appelle l'ouverture de la fenêtre d'ajout de données.
    /// </summary>
    public void OuvertureDuDialogAddData()
    {
        IServicesDocumentMxd service = FabriqueServices.Instance.CreerPersisteur<IServicesDocumentMxd>();
        service.OuvrirAddDataDialog();
    }
}

/// <summary>
/// Persisteur pour les objets de type mxDocument
/// </summary>
internal class ServicesDocumentMxd : IServicesDocumentMxd
{
     /// <summary>
     /// Ouverture de la fenêtre d'ajout de données.
     /// </summary>
     public void OuvrirAddDataDialog()
     {
          IAddDataDialog dialog = new AddDataDialog();
          bool ouvert = false;
          dialog.Document = (IMxDocument)this.applicationInterne.Document;
          dialog.Map = (this.applicationInterne.Document as IMxDocument).Maps.get_Item(0);
          dialog.Show(this.applicationInterne.hWnd, ouvert);
     }
}

The code above actually open's up the Add Data dialog, i need to figure out how to gain access to any .gdb to look for the layer inside them.

EDIT: Me

It actually works with the code above, as for the first edit, I have'nt tested it out at that moment wich i just did... Just wanted to share the progress... I will then proceed to close this question later today.