texts:
I am working on launching the GenerateNearTable using the ArcObjects Assembly method.  I have the advanced license (10.1), so the tool is available.  The code to launch geoprocessing is:

   Private Sub generateNearTable(ByRef FeatureLayer)

    'Grab the feature class for the GP tool
    Dim lcount As Long = 10
    Dim closest As Boolean = False

    'Setup the geoprocessor
    Dim gp As ESRI.ArcGIS.Geoprocessor.Geoprocessor = New ESRI.ArcGIS.Geoprocessor.Geoprocessor
    Dim neartable As GenerateNearTable = New GenerateNearTable()

    neartable.in_features = FeatureLayer
    neartable.near_features = FeatureLayer
    neartable.out_table = "C:\Users\jlaura\Desktop\test4" ' Hard code to test.
    neartable.closest = closest
    neartable.closest_count = lcount 'This should be a parameter

    Dim addtomap As Boolean = gp.AddOutputsToMap
    gp.SetEnvironmentValue("workspace", "C:\Users\jlaura\Desktop\test")
    gp.TemporaryMapLayers = False
    gp.AddOutputsToMap = False

    Try
        gp.Execute(neartable, Nothing)
        returnmessages(gp)
    Catch ex As Exception
        returnmessages(gp)
    End Try

End Sub

This is causing an error 99999.  The results are being returned to the results tab, so I am able to double click on the failed tool and run it successfully via ArcMap.  New to ArcObjects, so I'm definitely missing something obvious.

Do I need to check out a license like I would with ArcPy?
Am I parameterizing in a strange way?
Do I need to pass the dataset (IDataset) instead of the featureclass? #This failed when I tried it, but perhaps I missed something...

As an extension (that I can open as a new question if people think it is too disconnected) - to programatically set the output table, can I generate an empty, 'fieldless' standalone table and let the tool handle writing.

Updated Code:

Private Sub generateNearTable(ByVal distlayer, ByVal knn, ByVal DistanceTableOut)

    Dim FeatureLayer = GetFLayerByName(distlayer)
    Dim NearFeatureLayer = GetFLayerByName(distlayer) 'Fails if I use this.

    'Setup Params
    Dim lcount As Long = knn
    Dim closest As Boolean = False

    'Setup the geoprocessor
    Dim gp As New ESRI.ArcGIS.Geoprocessor.Geoprocessor
    Dim genneartable As GenerateNearTable = New GenerateNearTable()

    'Setup the Gp Environment
    Dim addtomap As Boolean = gp.AddOutputsToMap
    gp.TemporaryMapLayers = False
    gp.AddOutputsToMap = True 'Add the derived table to the map doc.
    gp.OverwriteOutput = True

    'Set Params
    genneartable.in_features = FeatureLayer
    genneartable.near_features = "C:\Users\jlaura\Desktop\Zunil\Zunil_secondaries_proj.shp"
    genneartable.out_table = DistanceTableOut
    genneartable.closest = closest
    genneartable.closest_count = lcount 'This should be a parameter

    'Execute
    Try
        gp.Execute(genneartable, Nothing)
        If gp.MessageCount > 0 Then
            For count As Integer = 0 To gp.MessageCount - 1
                MsgBox(gp.GetMessage(count), MsgBoxStyle.OkCancel, "Working")
            Next
        End If
    Catch ex As Exception

    End Try

End Sub

EDIT:

I am doing something wrong with the feature layer inputs.  By hard coding the path to the shapefiles I am able to successfully run the tool.  Back to the drawing board, re: featurelayer access.

EDIT2: The tools runs properly if the featurelayer is passed for the in_features and hard coded for the near_features.  In this case should I be reaccessing the feature layer from the input?


texts: Update polygon layer with summarized areas from other polygon layer

texts:
I use QGIS 3.2.3. I have two polygon layers "municipalities" and "landuse".

I want to update municipalities with the summarized area off landusepolygons within every municipality. landuse polygons do, in some cases, cross the borders of municipalities, in those cases I want only the intersection area in the calculation.

I can't find a way to do it with the build ind tools in QGIS. I have installed the plugin "refFunctions", but I can't figure out how to set it up.

A virtual layer with some sql would do it, I think. I can calculate the total area, but how do I get an summarized area for the landusepolygons in every municipality?


texts: Accuracy assessment using Google Earth Engine?

texts: I'm using Google Earth Engine for Land Cover change detection, I would like to know whether we can generate Error Matrix (Accuracy Assessment) for a classified image in Earth Engine or not?


texts: Get domain descriptions for multiple fields using SQL

texts:
How can I get domain descriptions for multiple fields using SQL?

I have successfully used the Oracle SQL in this ESRI document: Resolving domain codes to description values using SQL. But how do I do this for multiple fields?

SELECT
    EVENT_ID,
    ROAD_ID,
    Event_Status_Description
FROM
    TEMP_OWNER.ROADS
    LEFT OUTER JOIN
        (SELECT
           EXTRACTVALUE(CodedValues.COLUMN_VALUE, 'CodedValue/Code') AS Event_Status_Code,
           EXTRACTVALUE(CodedValues.COLUMN_VALUE,'CodedValue/Name') AS Event_Status_Description
        FROM
           SDE.GDB_ITEMS_VW items
           INNER JOIN SDE.GDB_ITEMTYPES itemtypes
           ON items.Type = itemtypes.UUID,
           TABLE(XMLSEQUENCE(XMLType(Definition).Extract ('/GPCodedValueDomain2/CodedValues/CodedValue')))     CodedValues
        WHERE itemtypes.Name = 'Coded Value Domain'
           AND items.Name = 'EVENT_STATUS') CodedValues
           ON TEMP_OWNER.ROADS.EVENT_STATUS = CodedValues.Event_Status_Code


texts: OSM download with overpass fails "Network request time out"

texts:
When trying to download OSM data in QGIS 2.8 I am using a layer to provide the extend for the requested data. The extend basically covers all the countries that share coast line of the Baltic sea (Russia only a small bit).

Now when I start the download process the process always fails with:
 Network request http://overpass-api.de/api/interpreter?data=(node(47.2701,5.86667,70.0924,41.3899);<;);out; time out

My questions:

is the extend to big?
is the amount of data requested to big?
is there an error in the network request?
is the function faulty in QGIS 2.8?
is there another way to extract OSM data, based on a given spatial extent?


texts:
I have installed overpass api on debian following install wiki page every thing is ok before starting the web api, when I start dispatcher daemon, It starts with no error and apache2 is started correctly but /api/interpreter is not found!404 Not Found, according to wiki page I tried :

sudo a2enmod cgi

but it says:

"Your MPM seems to be threaded. Selecting cgid instead of cgi."

Can anybody help me please?


texts: Set position of attribute table extended to next page

texts:
I am creating a script to export an attribute table in pdf format with PyQGIS 3. The number of features is different regarding the farmer.
pdf_table = QgsLayoutItemAttributeTable.create(layout)
pdf_table.setVectorLayer(pdf)
pdf_table.setMaximumNumberOfFeatures(250) # maximum number of features expected
layout.addMultiFrame(pdf_table)
frame2 = QgsLayoutFrame(layout, pdf_table)
frame2.attemptResize(QgsLayoutSize(198.200, 224),True)
frame2.attemptMove(QgsLayoutPoint(9.334, 31.071, QgsUnitTypes.LayoutMillimeters))
pdf_table.addFrame(frame2)
pdf_table.setResizeMode(1) # activate the extension of attribute table to the next page
#if the number of lines exceed the heigth defined

I would like to set the position of the attribute table extended in the next pages. Here is an example with an attribute table extended on one other page, but it can goes until 5 pages.

EDIT
An attribute table frame is created for each extension of the attribute table, so I need to get the reference of the extended ones and apply new parameters.

So I think it should be something like that :
for f in range(1, pdf_table.frameCount()): #start from 2nd page to the total number of sub-frames
    nextFrame = pdf_table.frame(f) #refers to the frame at index f
    nextFrame.attemptResize(QgsLayoutSize(198.200, 250),True)
    nextFrame.attemptMove(QgsLayoutPoint(9, 15), page=f)


texts: Map contents don't print when using CSS media print queries in a Bootstrap modal

texts:
My goal is to show an ArcGIS Server JavaScript v4 map in a Bootstrap modal popup (displayed over a page containing other content), and to print only the contents of the Bootstrap modal popup.

I'm using the approach at https://stackoverflow.com/a/12181825/906814 to successfully display the Bootstrap modal and print only its contents.

The issue is that while the ArcGIS map frame, attribution text and buttons print, the map's contents do not. Instead there is a blank section where the map contents should be:

What do I need to change in either the CSS or JavaScript to ensure that the map's contents print?

@media print {
  body * {
    visibility:hidden;
  }
  #printSection, #printSection * {
    visibility:visible;
  }
  #printSection {
    position:absolute;
    left:0;
    top:0;
  }
}

function printElement(elem, append, delimiter) {
   // from https://stackoverflow.com/a/12181825/906814
    var domClone = elem.cloneNode(true);

    var $printSection = document.getElementById("printSection");

    if (!$printSection) {
        $printSection = document.createElement("div");
        $printSection.id = "printSection";
        document.body.appendChild($printSection);
    }

    if (append !== true) {
        $printSection.innerHTML = "";
    }

    else if (append === true) {
        if (typeof (delimiter) === "string") {
            $printSection.innerHTML += delimiter;
        }
        else if (typeof (delimiter) === "object") {
            $printSection.appendChild(delimiter);
        }
    }

    $printSection.appendChild(domClone);
    }

Working demo at https://jsfiddle.net/slead/39Lrpbu8/3/ - press the Open The Modal button to display the map, then press the Print button in the modal's footer to display the print dialog.


texts:
I am attempting to interpolate sediments across an area with data from 3 sources using ArcGIS 10.3.

All three sources divide sediments into percentages of mud, sand or gravel.
Two sources come from actual sediment cores where grain size was analyzed in the lab.
The third comes from annotated photographs where a person estimates the % cover of each sediment type while viewing a photo of the seafloor.  The annotated photographs are understandably much less precise than the grain size analysis.

I have more data points from the annotated photos than the grain size analysis.

All percentages were converted into decimals before using geostatistical analyst to interpolate using ordinary kriging.
I will be repeating these steps for each sediment type (mud, sand, gravel).

Is it more appropriate to combine all data points into one data set and krig together, or is it better to use co-kriging and add each data-set as a separate layer/variable?

It seems like a gray area since the data are all either percent mud, sand or gravel but since some are collected with different methods, I could see using them as different variables.


texts: ArcMap is exporting to extremely large pdf files

texts:
My organization has one user who recently has noticed very large pdf files being exported from her mxds. When I recreate and do the export, I am getting a 4mb file (normal). When she exports the file we are seeing a 2.4 GB file. I changed some settings and was able to get it down to 1GB, but that wasn't satisfactory (obviously).

Things we have recently tried:

-Deleted the normal.mxt (reduced pdf size to 100mb)

-Turned off the setting that saves the spatial information in the pdf

-DPI is at 300

It isn't an overly large mxd. Normal vector data with some labeling. The pdf size is 42" x 60".

Any ideas on what to try next? She is on an ArcEditor 9.3.1 license.


texts: A small map take so much place and time?

texts:
I am using TileMill and try to generate an MBTile from Lyon in France.
Doing so takes me +100GB and an incalculable time (25mb for a really small portion of Lyon take me 1h)

If that is normal, does that mean services like Google Map store thousands of terabytes of images ?


texts: how to integrate geomajas i.e. clientside gwt base framework and serverside

texts: I am new to geomajas opensourse project.I have successfully build client-GWT maven project and server-side maven project.Now i want to integrate client side with server side project so that i can use server side geomajas as back end for client-GWT .I am not getting specific documents or commands to integrate client-side and server side. So anyone having knowledge of geomajas integration please provide me specific steps or commands ,so that i can successfully integrate client side and server side geomajas.


texts: What is the best method calculating travel distance between destination and point layer in QGIS when using OSGB 1936 / British National Grid coordinate reference system?


texts: Changing Points into Hexagon-Polygons using QGIS

texts:
I have only just started using QGIS and version 3.18, I need to convert Vector points into Polygons for a client and there was a plugin that did this in version 2 but I can't find anywhere in the questions that shows me how to do it in 3.18.
The answer below created the polygon for me but it only gives me three types of areas, rectangles, circles and diamonds, the client wants the shape to be a hexagon based on the width field in the data. The data field has a metre amount as in 5 for 5 metres.


texts: QGIS Redirecting WMS Request to localhost

texts:
Whenever I add a WMS layer, from a server on the local intranet, the WMS request gets sent to localhost instead of the actual server I am referencing.

I am enclosing screenshots to make sure my problem is clear.

Here's how I add a new WMS connection:

When I hit connect, the layers from the server are actually listed as shown below, showing that atleast the GetCapabilities request is going through fine:

However, the map doesn't show up and I get a 'Connection Refused' label under the map layer:

When I look in the logs, I see that the error message references 'localhost' instead of the server I set:

Clearly, something seems off. If I just copy the request from the log window and replace the localhost with my server IP into Firefox, then the request goes through fine and I see the map in Firefox.

Stranger still, is that if I look at the layer information, the server URL is referenced correctly:

This seems to be happening only in Windows QGIS builds. On Linux, the map shows up as expected.

Here is the QGIS version information:

I am posting here to check if I am doing something wrong. I can't seem to find a QGIS bug that describes this issue so I want to rule out that I am getting something wrong.


texts: Multiple rename shapefiles and feature classes in python

texts:
I have multiple shapefiles, in different sources (D:\GIS\Folder A", "D:\GIS_Temp\Name". And some of them have space in their name (File X, File Y, etc).

What should be the code that could rename all shapefiles in those folders, and change space with underscore, so the name of shapefiles would be File_X, File_Y, and so on...

When I try with this code:

import os, sys, arcpy
InWork =  ["D:\\GIS_Temp\Folder A", "D:\\GIS_Temp\\Folder B", "D:\\GIS_Temp\\Folder C"]
arcpy.env.workspace = InWork
for thisFC in arcpy.ListFeatureClasses():
newName = thisFC.replace(" ","_")
arcpy.Rename_management(thisFC,newName)

I put the folders in InWork, and I still get an error.

  Runtime error
  Traceback (most recent call last):
    File "", line 3, in
    File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\geoprocessing_base.py", line 515, in set_
      self[env] = val
    File "c:\program files\arcgis\desktop10.1\arcpy\arcpy\geoprocessing_base.py", line 567, in setitem
      ret_ = setattr(self._gp, item, value)
  RuntimeError: Object: Error in accessing environment

What am I doing wrong?


texts: Attempting to Draw Electricity Grid - Need Maximum Segment Length and Precise GPS Coordinates

texts:
I am attempting to draw an electricity grid (poles and cables) using QGIS. I need to ensure that each pole is a maximum of 50m apart, and then I need to measure the total length of the grid and identify the precise location of all poles. See my screenshot below for what it looks like - the yellow lines are electricity cable, the red dots are poles, and blue dots are houses.
I'm using the Advanced Digitizing tool to limit my line segments to be <=50m and Vector > Geometry > Extract Vertices to generate the pole locations. However, when I look at the attribute table for my poles shapefile, there are no GPS coordinates attached to my red dots (poles). This information is critical. How can I get GPS coordinates (ideally lat/long in decimal degrees) for the poles?
Is there a better way for me to do this? I need to draw a utility grid where no line segment exceeds 50m, and I need precise GPS coordinates of each pole.


texts: Methods in Slope, Aspect SAGA calculations

texts: Using this algorithm in QGIS I can choose between 7 methods, the default one being by Zevenbergen & Thorne 1987. What's the difference between them and can I just use the default one without worrying about it?


texts: Creating a virtual layer based on a virtual field for dynamic selection in QGIS 3.4?

texts:
Working on a PostGIS database, I would like to send to the print composer the selected features only. I'v tried several methods including working with is_selected() expression but the only way that actually works is constructing a filter on the table based on a list of values from attribute data which is quite tiresome.

I'v been trying to work with a virtual layer which could be dynamically updated with selected features from my source layer and that I could send to the print composer. I first set a virtual field on my source layer with the expression is_selected which returns 1 when the feature is selected. Then, to build the virtual layer, I first embed my source layer and then use a query on it :

SELECT * from my_layer where "selected" = '1'

"Selected" being my virtual field which is updated each time I select a new set of features.

However, even if the virtual layer toolbox is telling me there are no errors, the created virtual layer stays empty even if I have some selected features and my virtual field is returning '1' value. There are no values shown in the attribute table and no geometry shown on map.

Maybe working with a virtual layer is not the way to do it and there is a simpler way, but I'v not found it yet.


texts: FME Workbench Reader select tables to import using wildcards

texts:
In FME Workbench for various database format readers one can coose which tables to import:

Because all in all I have thousands of tables, I would like to select tables that begin with certain characters and in addition contain certain characters, for the example shown in the screenshot, I would like to select all tables beginning with 'strasseninfrastruktur' and additionally containing '_ct'.

Is there a way to do this?


texts: How to do a "continous display" in QGIS like in TileMill?

texts: In TileMill the "Earth" was automatically displayed "continuously", which means there was no "end". I tried to find this option in QGIS now, but nor search not checking through the options, etc. brought any success


texts: Add a delimited layer in QGIS (.csv), but do not add all columns

texts:
I am using the odk app to collect data and a heading of the form is required.

This heading is imported in QGIS as column as well, when I add the data as delimited layer, but it does not contain relevant information. Is there a way to remove certain columns while adding the .csv-file (skiping empty fields does not work)?

Image 1: xls-form (for ODK)

Image 2: Attribute Table in QGIS

I can delete the columns manuelly afterwards. But would like to know if there is a way to do it automatically (so either in the xls-form or before importing data in the settings of "adding a delimited text layer").


texts: How to batch geocode country names to polygons?

texts:
I have .xlsx dataset file with column "Country Name" and "Country Code" which have text records. I would like to geocode (assign) it to adequate country polygon in batch mode. The expected result should be a polygon of the country with attributes as columns from .xlsx.

Carto.com geocoder has such function (country geocode), however it is strictly limited to couple hundreds records and I have more than 6000 of it.

I am able to use QGIS and ArcGIS.


texts: Casting problem at CSWClientDemo of geotoolkit

texts:
I tried  to test CSWClientDemo of geotoolkit but couldn't succeed.

Here is a casting problem:

  GetRecordsResponseType response = ((JAXBElement<GetRecordsResponseType>) um.unmarshal(is)).getValue();

The error I got is:

Exception in thread "main" java.lang.ClassCastException: org.geotoolkit.ows.xml.v100.ExceptionReport cannot be cast to javax.xml.bind.JAXBElement
    at test.CSWClientDemo.main(CSWClientDemo.java:71)

(meanwhile, I imported org.apache.sis.xml.MarshallerPool instead of org.geotoolkit.xml.MarshallerPool.)


texts:
I followed the tutorial from http://schoolofdata.org/2013/04/27/creating-a-map-using-qgis  The France map is from NaturalEarthData and tour de france data is from another source.  They used different CRS (i.e., ED50 and WGS84 respectively).

I was wondering why the map and the data has different CRS?  Is this usual in GIS?


texts:
I am using OpenLayers 6 with GeoServer, and layers from my own creation with some shapefiles that were imported into a PostgreSQL database.
The idea is that the user clicks on the map, and the Points are saved on the database. I am trying to use WFS Transaction to achieve this, and following an example I created the following function:
function saveData(coordinate){
    var formatWFS = new ol.format.WFS();

    var formatGML = new ol.format.GML({
        featureNS: 'https://gsx.geolytix.net/geoserver/geolytix_wfs',
        featureType: 'wfs_geom',
        srsName: 'EPSG:3857'
    });

    var xs = new XMLSerializer();

    var iconFeature = new ol.Feature({
            geometry: new ol.geom.Point(coordinate),
            name: 'Null Island',
            population: 4000,
            rainfall: 500
        });
    var vectorSource = new ol.source.Vector({
            features: [iconFeature]
        });
    var dinamicPinLayer = new ol.layer.Vector({
            source: vectorSource
        });

    var node = formatWFS.writeTransaction(dinamicPinLayer, null, null, {
        featureNS: my_url,
        featureType: featureType,
        srsName: 'EPSG:3857'
    });

   //Also tried with:
  var node = formatWFS.writeTransaction(dinamicPinLayer, null, null, formatGML);

    var payload = xs.serializeToString(node);

    $.ajax('https://gsx.geolytix.net/geoserver/geolytix_wfs/ows', {
        service: 'WFS',
        type: 'POST',
        dataType: 'xml',
        processData: false,
        contentType: 'text/xml',
        data: payload
    }).done(function(response) {
        //sourceWFS.clear();
    });

}

When the user clicks on the map, the point is send to this function and the XML is sent to server. A 200 status is returned, but nothing is saved. When looking into the Request XML generated I only see the next Transaction:
<Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"/>

It should have an Insert tag with properties, but it only writes the transaction tag, so the GeoServer does not save the data. I tried different examples but all of them show this Transaction with no insert (nor update or delete) tag.
What am I missing?


texts: Select by expression functions in QGIS 2.8.2 (Windows)?

texts: Please forgive me if this is obvious, but having opened a newly installed version of QGIS on Windows, I can't seem to see the list of functions normally present (as shown in the Mac image) in the 'Select by expression' box. Anyone know how to get these back?


texts: FME: how to extract the coordinate of the North West corner

texts: In an FME translation from a polygon geodatabase feature class to a DGN I am trying to place a label in the top left corner of a box. I have tried using CoordinateExtractor to find the coordinates of a corner (then 2DPointReplacer and TexAdder to add text to this coordinate) but this works on the number of the vertex, which unfortunately isn't constant as not all of my boxes are quadrilateral. What I need to find is a way to find the upper left vertex. Does anyone know how to do this?


texts: Locate features along route - differentiate between real intersection and just touching bufferd zone Arcgis 10.4

texts:
Image: Case1 and Case 2

I have a polyline (also as route) representing a pipeline (green), a polyline representing the axis of a power line (black) and a polygon representing the allowable distance between powerline and pipeline (brown).

I use locate features along the route to identify the intersection of pipeline and powerline polygon to see where the pipeline is too close to the powerline. Result is a table with from x m to y m (FMEAS TMEAS) pipeline is in the polygon.
BUT:
If the pipeline is parallel to the power line the minimum distance (polygon) is compulsory.
If the pipeline intersects the powerline it is no problem.

How can I differentiate automatically between case 1 (not a real intersection with the powerline) and case 2 (a real intersection with the powerline)?


texts:
Currently trying to use gdal_merge within qgis to stack three raster bands, the command it's showing is like this:

  gdal_merge.bat -separate -of GTiff -o out.tif "LC80190302014126LGN00_B4.TIF" "LC80190302014126LGN00_B3.TIF" "LC80190302014126LGN00_B2.TIF"

however it keeps giving me ERROR 1: TIFFOpen:out.tif: Permission denied

I realize this should be a fairly simple procedure, is there something obvious I'm missing here?


texts: Is it possible to query via bounding box in the Google Maps Static API?

texts:
We have a Qt program written in C++ to request wms tiles and display them. This works fine.

Now we want to display imaginery from Google Maps and due to some restrictions in our program, we can't use javascript and we have to directly manage the url. For example, in the wms case:

www.idee.es/wms/PNOA-MR/PNOA-MR?SERVICE=WMS&VERSION=1.3.0&FORMAT=image/jpeg&CRS=epsg:32629&LAYERS=PNOA-MR&STYLES=&REQUEST=GetMap&EXCEPTIONS=INIMAGE&TRANSPARENT=TRUE&WIDTH=256&HEIGHT=256&BBOX=621921.34,4794536.07,622951.58,4801396.70

So I think the right choice for Google Maps would be the static API.

Now, reading it I didn't find anything regarding querying by bounding box. Only querying via a center and a zoom level.

Is that the only way to retrieve an image? Wouldn't it possible to query with a bounding box without specifying a center and the zoom level and retrieve an image representing that extension? Just like with the wms case before.

NOTE: I would like and url example if possible. We can't use javascript for this.


texts: LAS dataset to raster: 'no spatial reference exists'?

texts:
I am trying to make a DEM from my LAS dataset and I am using LAS Dataset to Raster tool in Conversion Tools. But I am getting the following error:

  No Spatial Reference exists.

  The table was not found.[VAT_Suwannee_Las4]

  ERROR 999999

The complete syntax submitted to the geoprocessing framework was:

arcpy.LasDatasetToRaster_conversion(in_las_dataset="C:/Users/Arsalan/Documents/ArcGIS/Suwannee.lasd", out_raster="C:/Users/Arsalan/Documents/ArcGIS/Default.gdb/Suwannee_Las3", value_field="ELEVATION", interpolation_type="BINNING AVERAGE LINEAR", data_type="INT", sampling_type="CELLSIZE", sampling_value="1", z_factor="1")

This is the screenshot showing the error:

Do you have any idea why this is happening?


texts: Failure to connect QGIS to ArcGIS services using WMS connector

texts:
I`m trying to add a layer to QGIS connecting to a WMS server but I got this error:

  "The server you are trying to connect to is not a WMS server. Please check the URL"

I tried two different URLs but I got the same message. Here are the URLs:

http://www.snirh.gov.br/arcgis/rest/services/SNIRH2016/Precipitacao_isoietas/MapServer/1

and

https://services.arcgisonline.com/arcgis/rest/services/Elevation/World_Hillshade/MapServer

What I'm doing wrong? Sorry for the newbie question.

Step 1:

Step 2:

Here is the detailed message:

"Em vez da string de caracteres de recursos que era esperada, foi recebida a seguinte resposta: (in portuguese)"

translated:

"Instead of the resource string that was expected, the following response was received"


texts:
I am using QGIS 2.4.0 and sadly it seems that QGIS imported csv files are somewhat bugged. I cannot use imported csv files the same way as shapefiles. It seems that sometimes it is easier to import the file and then just save the csv file inside QGIS as a shapefile and import it again to make it work.
The reasons I am saying this are the following:

I cannot save the style files associates with csv files.
In properties tab type and typename are always double and the lenght and precision are always 0
Calculations are not possible since some columns are not seen as containing Real or Natural numbers.

Once saved as shapefile I can suddenly do everything I need. So am I missing something or is there no other way to not to have this intermediary saving step?


texts:
The ArcGIS Javascript API has a nice refresh interface for ArcGISDynamicMapServiceLayer:

    var dynamicMapServiceLayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Demographics/ESRI_Census_USA/MapServer");
    map.addLayer(dynamicMapServiceLayer);
    dynamicMapServiceLayer.setDisableClientCaching(true);

    if (dynamicMapServiceLayer.disableClientCaching) {
        aspect.after(dynamicMapServiceLayer, "onUpdateEnd", function () {
            $.delay(15000, "UpdateDynamicLayer", function () {
                dynamicMapServiceLayer.refresh();
            });
        });
    }

But this is not generalized across the other layer types.  Anyone have a solution for refreshing layers of any type?


texts:
I use a Slak 14 box and tried to install QGIS on it. the compile process goes on for 30 minutes or so (dual Xeon 3.2 ghz) and then I get the errors below. Says the compiler doesn't work. Have had no other problem with this compiler. Anybody able to point me in a direction? I'm kinda not sure which way to pursue this.The one thing strikes me is that it cmake isn't acting well.

...many minutes of output.

  qgis-2.4.0/tests/testdata/zonalstatistics/polys.cpg
  qgis-2.4.0/tests/testdata/zonalstatistics/polys.dbf
  qgis-2.4.0/tests/testdata/zonalstatistics/polys.prj
  qgis-2.4.0/tests/testdata/zonalstatistics/polys.qpj
  qgis-2.4.0/tests/testdata/zonalstatistics/polys.shp
  qgis-2.4.0/tests/testdata/zonalstatistics/polys.shx
  qgis-2.4.0/tools/
  qgis-2.4.0/tools/qgis-qt3to4.pl
  qgis-2.4.0/tools/qgis_config/
  qgis-2.4.0/tools/qgis_config/qgis_config.cpp
  qgis-2.4.0/tools/qgisdetect.m4
  -- The C compiler identification is unknown
  -- The CXX compiler identification is unknown
  -- Check for working C compiler: /usr/bin/cc
  -- Check for working C compiler: /usr/bin/cc -- broken
  CMake Error at /usr/share/cmake-2.8/Modules/CMakeTestCCompiler.cmake:61 (message):
    The C compiler "/usr/bin/cc" is not able to compile a simple test program.

It fails with the following output:

  Change Dir: /tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp

  Run Build Command:/usr/bin/gmake "cmTryCompileExec2321741014/fast"

  /usr/bin/gmake -f CMakeFiles/cmTryCompileExec2321741014.dir/build.make
    CMakeFiles/cmTryCompileExec2321741014.dir/build

  gmake[1]: Entering directory
    `/tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp'

  /usr/bin/cmake -E cmake_progress_report
    /tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp/CMakeFiles 1

  Building C object
    CMakeFiles/cmTryCompileExec2321741014.dir/testCCompiler.c.o

  /usr/bin/cc -O2 -march=i486 -mtune=i686 -o
    CMakeFiles/cmTryCompileExec2321741014.dir/testCCompiler.c.o -c
    /tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp/testCCompiler.c

  /tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp/testCCompiler.c:1:0: error:
    CPU you selected does not support x86-64 instruction set

  #ifdef __cplusplus
     ^

  gmake[1]: Leaving directory `/tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeTmp'

  gmake[1]: * [CMakeFiles/cmTryCompileExec2321741014.dir/testCCompiler.c.o]
    Error 1

  gmake: * [cmTryCompileExec2321741014/fast] Error 2

  CMake will not be able to correctly generate this project.
  Call Stack (most recent call first):
    CMakeLists.txt:7 (PROJECT)

  -- Configuring incomplete, errors occurred!
  See also "/tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeOutput.log".
  See also "/tmp/jah/qgis-2.4.0/build/CMakeFiles/CMakeError.log".
  Cannot install /tmp/qgis-2.4.0-i486-1-jah.txz:  file not found

  Installation complete

  bash-4.2#


texts: ST_Touches for more than just one geometry

texts:
I'd like to select some neighbor polygons with ST_Touches. Therefore I found the following statement:

SELECT b.name
FROM layera as a
JOIN layera as b
  ON ST_Touches((SELECT a.geom FROM layera as a WHERE a.id = 5124),b.geom)
GROUP BY b.name

With this I get all the neighbor polygons of the polygon with the id '5124'. But now I want to call this function with a list of arguments. See the example statment below:

SELECT b.name
FROM layera as a
JOIN layera as b
  ON ST_Touches((SELECT a.geom FROM layera as a WHERE a.id > 5000),b.geom)
GROUP BY b.name

(SELECT a.geom FROM layera as a WHERE a.id > 5000 returns lets say 100 results)

I know that this won't work because ST_Touches works for only two arguments. But this should show what I'd like to do. How is it possible to insert all the results of a SELECT-statement  into this ST_Touches function?

Thanks for any suggestions. =)

EDIT:

Ok, I just found out that I get the same result set with this statement:

SELECT layera.name
FROM layera
WHERE ST_Touches(layera.geom,(SELECT layera.geom FROM layera WHERE layera.nr = 5134))

Anyway, the problem isn't solved. =(


texts: Charting two FeatureCollections in Google Earth Engine?

texts:
I have 2 shapefiles representing different types of land-use. Each Shapefile contains 10 polygons. I would like to chart the mean NDVI value over a specific timespan for each class. How can I display 2 FeatureCollection in the same chart?

I tried .flatten(), but then the chart displays only 1 graph.

var L8 = ee.ImageCollection("LANDSAT/LC8_L1T_TOA"),
extensiv = ee.FeatureCollection("users/tobkugler/Engadin_Extensive_Merc_new"),
intensiv = ee.FeatureCollection("users/tobkugler/Engadin_Intensiv_Merc_new");

var classes = new ee.FeatureCollection([extensiv,intensiv]);
var classes_flatten = classes.flatten();

// Timespan for the ImageCollection

 var filtered = L8.filterDate('2014-01-01', '2016-12-31')
.filterBounds(extensiv);

// Function to add NDVI band

function addNDVI(image) {
var ndvi = image.normalizedDifference(['B5', 'B4']);
return image.addBands(ndvi);
var with_ndvi = filtered.map(addNDVI);

//Print NDVI Chart

 var chart = ui.Chart.image.doySeriesByRegion(with_ndvi,'nd',classes_flatten,
 ee.Reducer.mean(),30,ee.Reducer.mean(),'label');

 print(chart);


texts: Difference between NAD83 Indiana SPW in QGIS and ArcGIS for Desktop?

texts:
It seems ESRI and QGIS don't have the same definition for my most commonly used projection, NAD83 StatePlane Indiana West (US ft).  When I create data in an ESRI product the projection shows up in QGIS as:

  USER:100001 -  * Generated CRS (+proj=tmerc +lat_0=37.5
  +lon_0=-87.08333333333333 +k=0.9999666666666667 +x_0=900000 +y_0=250000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=us-ft +no_defs)

When I create data in QGIS the projection information shows up under the Custom folder as:

NAD83_Indiana_West_ftUS
Authority: Custom

Projection: Transverse_Mercator
false_easting: 2952750.0
false_northing: 820208.333
central_meridian: -87.08333333333333
scale_factor: 0.999966667
latitude_of_origin: 37.5
Linear Unit: Foot_US (0.3048006096012192)

Geographic Coordinate System: GCS_North_American_1983
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_North_American_1983
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

ESRI's normal entry is:

NAD_1983_StatePlane_Indiana_West_FIPS_1302_Feet
WKID: 2966 Authority: EPSG

Projection: Transverse_Mercator
False_Easting: 2952750.0
False_Northing: 820208.3333333333
Central_Meridian: -87.08333333333333
Scale_Factor: 0.9999666666666667
Latitude_Of_Origin: 37.5
Linear Unit: Foot_US (0.3048006096012192)

Geographic Coordinate System: GCS_North_American_1983
Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_North_American_1983
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

How big of a deal is this?  The only difference I can see is in the precision of the scale factors.  Is there another difference I missed?  Is it worth worrying about?


texts: How to programmatically create a bivariate choropleth in d3.js?

texts:
Question and Situation
I'm using d3.js to visualize some data on a map. The user clicks on one checkbox and all path elements get filled according to a color-preset, which all works fine. Now what we also have decided is that the user can chose up to 2 concurrent variables which have to get displayed in a bivariate choropleth. It works on my side but the color mix looks somewhat mashed and singular (see screenshots) - How do i get a correct mix of these colorscales?
Code
const colorize = []
if (this.selectedVariables.length > 0) {
  const vars = this.selectedVariables
  vars.forEach(colorData => {
    const { legend } = colorData
    const color_range = legend.values.map(color => Number(color.value))
    colorize.push(d3.scaleThreshold(color_range, d3[`scheme${legend.color_scheme}`][color_range.length]))
  })
}

// The variables from this.selectedVariables
[{
    "group": "Variables",
    "name": "crime_rate",
    "type": "checkbox",
    "label": "Crime Rate",
    "description": "
Info about Crime Rate
",
    "legend": {
        "title": "Crime Rates",
        "color_scheme": "YlOrRd",
        "values": [{"value": "3891"}, {"value": "5258"}, {"value": "7091"}, {"value": "9962"}, {"value": "12821"}, {"value": "18832"}, {"value": "520436"}]
    }
}, {
    "group": "Variables",
    "name": "rel_prop",
    "type": "checkbox",
    "label": "Christianity",
    "description": "
Christianity by Rel Prop
",
    "legend": {
        "title": "Christianity ",
        "color_scheme": "BuPu",
        "values": [{"value": "0.21"}, {"value": "0.53"}, {"value": "0.59"}, {"value": "0.64"}, {"value": "0.69"}, {"value": "0.75"}, {"value": "0.85"}]
    }
}]

// Fill all path elements
pathData
  .attr('d', (d) => pathGen(d.geometry))
  .attr('fill', d => {
    if(this.selectedVariables.length > 0 && d.properties['mapData']) {
      if (colorize.length < 2) {
        return colorize[0](d.properties['mapData'][this.selectedVariables[0].name])
      } else {
        return d3.scaleLinear().range(colorize.map((color, i) => color(this.selectedVariables[i].name)))(0.5)
      }
    } else {
      return '#444'
    }
  })

Screenshots


texts: Error code 500 at REST Endpoint of published geoprocessing service perhaps due to output parameter type?

texts:
Using ArcGIS 10.1 (Server and Desktop) I've been trying to publish a geoprocessing service based on a python script I wrote. Of course, it all works flawlessly from ArcGIS Desktop but I run into all sorts of problems when publishing.

I'm pretty sure the problem is related to the types of my output parameters or something to that effect. Here's the python script, minus comments and messages to GUI...

import arcpy
import json
import sys
import os

from arcpy import env

class Result(object): pass

out = Result();

arcpy.env.overwriteOutput = True

config = arcpy.GetParameterAsText(0)    # A string containing a path to a JSON config file
source = arcpy.GetParameter(1)          # A string containing the name of a layer
sWhere = arcpy.GetParameter(2)          # A String containig a where clause
distan = arcpy.GetParameter(3)          # A String containing a distance linear unit
target = arcpy.GetParameter(4)          # An Array containing multiple layer names

try:
    f = open(config)

    try:
    jsonCfg = json.load(f)

    finally:
    f.close()

except:
    pass

# jsonCfg['connection'] is a path to an SDE connection
fc = os.path.join(jsonCfg['connection'], source)
scratch = os.path.abspath(jsonCfg['scratch'])

arcpy.MakeFeatureLayer_management(fc,"source_lyr")
arcpy.SelectLayerByAttribute_management("source_lyr", "NEW_SELECTION", sWhere)
bufLyrNam = os.path.join(scratch, "source_buf")
arcpy.Buffer_analysis("source_lyr", bufLyrNam, distan)

buffer = arcpy.FeatureSet()
buffer.load(bufLyrNam)

nLayers = len(target)

selection = []

for i in range(0, nLayers):
    lyrNam = "target_lyr_" + str(i);
    tgtNam = os.path.join(jsonCfg['connection'], target[i])

    arcpy.MakeFeatureLayer_management(tgtNam, lyrNam)
    arcpy.SelectLayerByLocation_management(lyrNam, "INTERSECT", bufLyrNam)
    cpyNam = os.path.join(scratch, target[i].replace(".",""))
    arcpy.CopyFeatures_management(lyrNam, cpyNam)

    fs = arcpy.FeatureSet()
    fs.load(cpyNam)
    selection.append(fs)

out.selection = selection

arcpy.SetParameter(5, out)

At first, I wanted to output the "out" object but apparently, this is not allowed in geoprocessing services. I can still publish the service but the output is converted to an unusable string.

I then tried to output only the "buffer" variable (containing a  featureset) using an output parameter of type Any Value, Feature Class, Feature Layer... etc. and tried to output the selection variable as a Multiple Value Feature Class with no success.

In fact, something weird happens when I try to output the buffer. Once my service is published, I cannot access its REST endpoint, I always get a error code 500 page. This leaves me thinking that the arcpy.SetParameter must correspond to the output parameter type in order to even access the REST endpoint.

I've been working on this for a couple of days and I'm out of ideas which is why I created an account and asked my question. Can anyone help me figure out the problem here?


texts: On click GeoServer WMS layer features in ArcGIS JS API

texts:
I am consuming a GeoServer WMS layer in ArcGIS JavaScript API. I am successfully able to display layer. Now I am not able to get layer info on map click or layer click.

Below is code using to add GeoServer WMS layer, how can I get the layerinfo?

var layer1 = new WMSLayerInfo({
    'id': '1',
    name: 'xxx:yyy',
    title: 'nrega',
    queryable: true,
    });
var resourceInfo = {
    featureInfoFormat: "text/html",
    infoTemplate: infoTemplate,
    tolerance: 10,
    extent: new Extent(7971379.79404854, 1129044.28557475, 9134620.92140827, 2217652.68573662, new SpatialReference({
        "wkid": 4326
        })),
    getFeatureInfoURL: "https://adcd/xyz/ows/?service=WMS&version=1.3.0&info_format=application/json&xsl_template=infoTemplate.xsl",
    layerInfos: [layer1]
    };
var wmsLayer = new WMSLayer("https://adcd/xyz/wms/", {
                                resourceInfo: resourceInfo,
                                visibleLayers: ['xxx:yyy'],
                                visible: true
    });
wmsLayer.spatialReferences[0] = 3857;
wmsLayers = wmsLayer;
map.addLayers([wmsLayers]);


texts:
My last question here was " Install webpage Admin in Ubuntu GIS (Server) " but i´m trying to do something new while i´m waiting for the possible answers:

I have Windows 10 in PC (With QGIS 2.18.7 and PostgreSQL 4v1 installed) and i installed Virtual Box to run UBUNTU Desktop 17.04 as you can see:

My question is: What is the best spatial data server I can use in ubuntu to allow connection to Postgresql?

NB: With a friendly work environment.

I have read it:

a) https://geoind.wordpress.com/2014/02/17/r-e-rstudio-no-ubuntu/

b) http://www.oslandia.com/full-spatial-database-power-in-2-lines-en.html

c) https://wiki.osgeo.org/wiki/Wms

d) https://help.ubuntu.com/community/UbuntuGIS


texts:
I have made a Button addin for ArcMap in Visual Studio 2010, but I don't know how to deploy it.

I tried debugging it from ArcMap but nothing happens.


texts: ArcGIS feature layer built of SQL View not visible

texts:
I am building a custom web app using the ArcGIS JS 4.8 API with an ArcGIS Web service and I have run into an issue.

The root cause that I am trying to fix is to do with features not being filtered out when I perform queries, due to them having duplicate ObjectID values.

I believe this is to do with left-joining non-spatial data to the spatial data since this is a time-based app (in essence, if a road or area requires maintenance every month, this app will have one spatial segment with multiple non-spatial segments, which duplicates the spatial unique ID).

The current table in our production environment looks like this, and when I perform queries on our data I will always see 10001, and 10002, even if I filter it to only show where objectID=10003

| OBJECTID | MAINTENANCEID |  MAINTENANCENAME  | DATE |
|==========|===============|===================|======|
| 10001    | 1             | Filter Sand in pit|2017  |
| 10001    | 2             | Filter Sand in pit|2018  |
| 10002    | 3             | Filter Sand in pit|2019  |
| 10002    | 4             | Mow field grass   |2018  |
| 10003    | 5             | Mow field grass   |2019  |

To fix this, I have created an SQL view of my spatial and non-spatial data with a composite unique ID made of the spatial and non-spatial unique id's, which is published through ArcMap as a Web Service. The new column is formatted very similarly to the NEW_OBJECTID field and is not a stored field but a dynamically generated identity value.

| OBJECTID | MAINTENANCEID |  MAINTENANCENAME  | DATE | NEW_OBJECTID |
|==========|===============|===================|======|==============|
| 10001    | 1             | Filter Sand in pit|2017  | 1000101      |
| 10001    | 2             | Filter Sand in pit|2018  | 1000102      |
| 10002    | 3             | Filter Sand in pit|2019  | 1000203      |
| 10002    | 4             | Mow field grass   |2018  | 1000204      |
| 10003    | 5             | Mow field grass   |2019  | 1000305      |

Unfortunately, when I try to display this web service data on a Map and MapView, I can not see the works using a JS FeatureLayer.

Inspecting the Web Service through the ArcGIS/rest/services/ online viewer, I can see and query the data. This, however, uses an ArcGIS API 3.25 ArcGISDynamicMapServiceLayer instead of an API 4.8  FeatureLayer.

The front-end code I am using is as simple as I could make it reduce errors, and I have tested it with the current production web service, and it shows correctly.

<script type="text/javascript">
require([
    "esri/Map",
    "esri/views/MapView",
    "esri/layers/FeatureLayer",
    "dojo/domReady!"
  ],
  function(
    Map, MapView,
    FeatureLayer
  ) {

  var map = new Map({
      basemap: "streets"
    });
  var view = new MapView({
    map: map,
    container:"map",
  });

    var layer = new FeatureLayer({url: "<MY_URL>"});
    map.add(layer);
    layer.when((x) => {console.log(x)})
});
</script>

Any insights on how I could debug or fix this?


texts:
I'm looking to be able to update a value by its nearest polygon.

The below image shows polygons in green where the column roadclass = 'A' and shows polygons where column roadclass IS NULL in brown.

I'm manually selecting and updating these currently so would like to automate it if possible.

I'm keen to learn but I've not used python before in QGIS so please consider this in any response. Also if useful the layer is stored within a PostGIS database and I'm using QGIS 3.6.


texts: Is it Better to Use Google Maps or OpenStreetMap When Dealing With Data Points in both Urban and Rural Areas?

texts:
My application tracks cancer-related trends across multiple US States.
Most of the patient data in the application deals with people who live in poor, rural, or under-served regions.
We want to start taking data from some of our demographic-related reports and plotting them on maps (e.g., pointing out cancer trends / screening by region and county)

That all being said, is it better for us to go with an Enterprise-level license of Google Maps due to the way the data is collected (by professionals with fancy tech)? Or is it better for me to go with OpenStreetMap since it may contain mapping data for more of the extremely remote areas some of our patients live in?


texts: How to get user history from OS maps API?

texts:
Is there a way to get the history of a user as found on the OS maps

by using the OS maps API? https://api.openstreetmap.org/api/0.6/user/8566797


texts: Convert OSM into shapefile with GeoTools?

texts: Is it possible to convert the OSM file to shapefile with GeoTools in Java?


texts: Sentinelsat: Error trying to download from Norweigan mirror colhub.met.no

texts:
I am trying to download the Sentinel 2 images over an area of Western Norway. Since the ESA archive is very slow, I am trying to use the Norwegian mirror for Sentinel instead (https://colhub.met.no/#/home).

I manage to get an excel file of the images suitable using the following code (username & password taken out):

from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
import pandas as pd
import logging

logging.basicConfig(format='%(message)s', level='INFO')
api = SentinelAPI('username', 'password', 'https://colhub.met.no/#/home')
dates = '[2018-01-01T00:00:00.000Z TO 2019-07-01T00:00:00.000Z]'
orbit = '117'
footprint = geojson_to_wkt(read_geojson("P:\Ben\hordaland_outline_simp.json"))
cloud = '[0 TO 50]'

products = api.query(footprint, dates,platformname='Sentinel-2',cloudcoverpercentage = cloud, producttype = 'S2MSI1C')

# convert to Pandas DataFrame
products_df = api.to_dataframe(products)
##products_df_sorted = products_df.sort_values(['beginposition', 'orbitdirection'], ascending=[True, True])

file_name = r"G:\Data\Raw\Sentinel\Sentinel2\1_Raw_L2A\Sentinel2_images.xlsx"
# Create a Pandas Excel writer using XlsxWriter as the engine.
writer = pd.ExcelWriter(file_name, engine='xlsxwriter')
# Convert the dataframe to an XlsxWriter Excel object.
products_df.to_excel(writer, sheet_name='Sheet1')
# Close the Pandas Excel writer and output the Excel file.
writer.save()
print "Excel file with images to be processed has been written to " + file_name

Then I try to download using

api.download_all(products_df)

But I get the error:

Traceback (most recent call last):
  File "C:\Users\bro076\AppData\Roaming\Python\Python27\site-packages\sentinelsat\sentinel.py", line 602, in download_all
product_info = self.download(product_id, directory_path, checksum)
  File "C:\Users\bro076\AppData\Roaming\Python\Python27\site-packages\sentinelsat\sentinel.py", line 493, in download
product_info = self.get_product_odata(id)
  File "C:\Users\bro076\AppData\Roaming\Python\Python27\site-packages\sentinelsat\sentinel.py", line 428, in get_product_odata
values = _parse_odata_response(response.json()['d'])
  File "C:\Users\bro076\AppData\Roaming\Python\Python27\site-packages\sentinelsat\sentinel.py", line 1128, in _parse_odata_response
'Online': product['Online'],
KeyError: 'Online'

I was wondering if some of the parameters need to be tweaked when downloading from colhub.met.no, perhaps since the data is not archived in the long-term storage (Online/Offline) as is the ESA archive.


texts: When does iface.mapCanvas().layers() get updated to the contents of QgsMapLayerRegistry.instance().mapLayers()?

texts:
Since the topic might be a bit cryptic, I'll try to explain where my problem is.

I am currently developing a Qgis-plugin, where I am creating and adding some layers automatically, which I would like to activate afterwards. So far, so good, not a big deal. Hence I am doing this process several times, I created a method for iterating over the layers.

And now the problem is:

I create and add the layer:

        type = 'Polygon'
        crs = '?crs=' + self.project_crs
        shape_layer = QgsVectorLayer(type + crs, layer_name, 'memory')
        shape_layer.setProviderEncoding('System')

        # add the layer to the layer-legend
        QgsMapLayerRegistry.instance().addMapLayer(shape_layer)

But if I want to access the layer by iterating over the iface.mapCanvas().layers() and running:

            for layer in self.iface.mapCanvas().layers():

                if layer.name() == layer_name:
                    return layer

'None' is returned. Strange thing, cause the problem could be avoided, by either accessing the mapCanvas() prior to starting iteration:

            layers = self.iface.mapCanvas().layers()
            for layer in self.iface.mapCanvas().layers():

                if layer.name() == layer_name:
                    return layer

Or, by iterating over either self.iface.legendInterface().layers() or QgsMapLayerRegistry.instance().mapLayers(), which both are "updated" instantly.

Could anyone explain the mechanisms behind the addMapLayer()-function and in general regarding the memory-mgmt behind Qgis? Or does anyone know where to read something about it?


texts: OGR DeleteDataSource not working properly in Python

texts:
I have a small Python script which is supposed to delete a shapefile in case the file is empty.

import sys
from osgeo import ogr

drv = ogr.GetDriverByName("ESRI Shapefile")
ds = drv.Open(inData)
if ds is None:
    print "Could not open file."
    sys.exit(1)

lyr = ds.GetLayer()

if lyr.GetFeatureCount() == 0:
    drv.DeleteDataSource(inData)

everything seems to work fine. The shapefile consists of

shapefile.dbf
shapefile.prj
shapefile.shp
shapefile.shx

whereas only the *.prj and *.shx get deleted. the *.dbf and *.shp are untouched. Is there anything I am missing?


texts: Raster Calculator QGIS 3.0 , cannot specify output data type

texts:
I am trying to generate a raster layer from the subtraction of a DSM raster layer and a DTM raster layer.
Both my DSM and DTM are in float64 data type. When I run the raster calculator to do the subtraction, I could not specify the data type of the output and the default output that I got was in float32 after subtraction. Is there an option and button to specify the output data type?

I am using QGIS 3.0.


texts: Copying Shape_Length and Shape_Area data (replace geometry tool) to other rows, or suitable alternatives using ArcMap?

texts:
I am working on mapping extents for historical construction plans using ArcMap 10.6. I have sets of these blueprints that I must draw out dozens of similar polygons using the replace geometry tool, but on the larger and more complex polygons, I find myself making errors that slow my workflow down.

Is there is not a way for me to copy the polygon data, so that I only need to make one polygon that can be moved to other rows, what would be an alternative other than just continuing to draw out each sheet?


texts: Select features within a country or bounding box effectively using PostGIS

texts:
I imported the planet OSM file into a PostgreSQL 9.2/PostGIS 1.5 server using osm2pgsql. I now want to select certain map features which lie within the bounding box formed by a countries borders, e.g. all rivers and lakes within Germany.

I've come up with a working query, but it's horribly slow:

osm=# EXPLAIN ANALYZE
WITH     bbox
AS        (SELECT ST_Envelope(ST_Union(way)) AS way
        FROM    planet_osm_polygon
        WHERE     boundary = 'administrative' AND admin_level = '2' AND name = 'Deutschland')
SELECT    *
FROM    (
            (SELECT     way, tags, width, water, waterway
            FROM        planet_osm_line
            WHERE waterway='river')
        UNION
            (SELECT     way, tags, width, water, waterway
            FROM        planet_osm_polygon
            WHERE        water IS NOT NULL)
        ) AS q, bbox
WHERE    ST_Within(q.way, bbox.way);
                                                                                QUERY PLAN
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 Nested Loop  (cost=1558037.52..1676764.41 rows=1 width=192) (actual time=305373.041..311176.612 rows=20435 loops=1)
   Join Filter: ((planet_osm_line.way && bbox.way) AND _st_within(planet_osm_line.way, bbox.way))
   Rows Removed by Join Filter: 851602
   CTE bbox
     ->  Aggregate  (cost=18.27..18.28 rows=1 width=1005) (actual time=374.222..374.222 rows=1 loops=1)
           ->  Index Scan using planet_osm_polygon_boundary_admin_level_name_idx on planet_osm_polygon  (cost=0.00..18.26 rows=1 width=1005) (actual time=27.989..34.350 rows=8 loops=1)
                 Index Cond: ((boundary = 'administrative'::text) AND (admin_level = '2'::text) AND (name = 'Deutschland'::text))
   ->  CTE Scan on bbox  (cost=0.00..0.02 rows=1 width=32) (actual time=374.227..374.228 rows=1 loops=1)
   ->  Unique  (cost=1558019.24..1564213.69 rows=412963 width=1131) (actual time=299869.930..309268.592 rows=872037 loops=1)
         ->  Sort  (cost=1558019.24..1559051.65 rows=412963 width=1131) (actual time=299869.929..300244.843 rows=874853 loops=1)
               Sort Key: planet_osm_line.way, planet_osm_line.tags, planet_osm_line.width, planet_osm_line.water, planet_osm_line.waterway
               Sort Method: external sort  Disk: 693392kB
               ->  Append  (cost=0.00..1314829.27 rows=412963 width=1131) (actual time=64.278..204999.538 rows=874853 loops=1)
                     ->  Index Scan using planet_osm_line_waterway_idx on planet_osm_line  (cost=0.00..440480.43 rows=132130 width=1245) (actual time=64.278..142495.404 rows=595834 loops=1)
                           Index Cond: (waterway = 'river'::text)
                     ->  Bitmap Heap Scan on planet_osm_polygon  (cost=5271.85..870219.21 rows=280833 width=1078) (actual time=186.014..62347.033 rows=279019 loops=1)
                           Recheck Cond: (water IS NOT NULL)
                           ->  Bitmap Index Scan on planet_osm_polygon_water_idx  (cost=0.00..5201.64 rows=280833 width=0) (actual time=123.545..123.545 rows=279019 loops=1)
                                 Index Cond: (water IS NOT NULL)
 Total runtime: 311350.392 ms

I've already touched some PostgreSQL settings:

shared_buffers = 2GB
temp_buffers = 16MB
work_mem = 32MB

I tried to get rid of the external disk sort with a high value for work_mem and set work_mem = 768MB, but the disk sort is still in the query plan.

Is there a way to speed things up?


texts: How to check if GeoJSON feature is in rectangular shape and find corner points?

texts:
I have a collection of GeoJSON features as Polygon & MultiPolygon which are saved in MongoDB. Many among them are in rectangular or square shape, while others are odd shaped. How can I find all the four corner (top left, top right, bottom left and bottom right) points if that feature is in rectangular or square shape?

I first tried to filter those features which has only five co-ordinates, so that those points will be corner co-ordinates essentially. But some of them are having more than five co-ordinates but as shape they are rectangle or square. Check below given feature examples. First is odd shaped feature, second is rectangle with redundant points and third is rectangle with 4+1 points.

I checked TurfJS, but didn't find any helpful method. BBox gives me non-rotated rectangle, which is not useful because it can also give bbox for non-rectangular shape also.

{
  "type": "FeatureCollection",
  "features": [
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                -102.61383274599996,
                32.37585257400008
              ],
              [
                -102.61764297999997,
                32.375000929000066
              ],
              [
                -102.62239479299996,
                32.388750307000066
              ],
              [
                -102.60775890199994,
                32.392205895000075
              ],
              [
                -102.60566771099997,
                32.392699906000075
              ],
              [
                -102.60510962499995,
                32.39093294600008
              ],
              [
                -102.61772064899998,
                32.387951290000046
              ],
              [
                -102.61383274599996,
                32.37585257400008
              ]
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                -102.59602204999999,
                32.36491866700004
              ],
              [
                -102.59952111499996,
                32.37500140000003
              ],
              [
                -102.60087906199999,
                32.378914276000046
              ],
              [
                -102.58443896599994,
                32.38292495700006
              ],
              [
                -102.58184811899997,
                32.37500015000006
              ],
              [
                -102.57980711899995,
                32.36875659300006
              ],
              [
                -102.59602204999999,
                32.36491866700004
              ]
            ]
          ]
        ]
      }
    },
    {
      "type": "Feature",
      "properties": {},
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                -102.62506373399998,
                32.34304188700003
              ],
              [
                -102.64139887699997,
                32.33901489700003
              ],
              [
                -102.64604796299994,
                32.352882700000066
              ],
              [
                -102.62962431499994,
                32.35679300500004
              ],
              [
                -102.62506373399998,
                32.34304188700003
              ]
            ]
          ]
        ]
      }
    }
  ]
}

https://bl.ocks.org/Xyroid/raw/800dafa55111ce139b8158c62f858c98/


texts:
We are planning to implement a routing application to monitor marine traffic.

The issue is that we either need a data set that has international sea traffic routes data or a clear definition of it.

In both cases, we couldn't find a good resource. We found few resources such as geoweb.dnv.org/data and http://www.arcgis.com/home/item.html?id=3d82f13e484c4d3ca24345b796b76400 but it's limited for few areas, while we need to cover the whole globe.

Any idea?

to explain what I'm looking for, please take a look here:
http://www.imo.org/OurWork/Safety/Navigation/Pages/ShipsRouteing.aspx

I found this dataset http://geocommons.com/overlays/25 but it's not quite accurate. I don't think it's based on the standard traffic routing system.


texts: Can't get Tiger Geocoder to work in PostGIS

texts:
I've installed the PostGIS Tiger Geocoder according to these instructions and loaded in data for the state of New York.

Unfortunately, when I try to geocode something...

ny_tiger=# SELECT g.rating, ST_X(g.geomout) As lon, ST_Y(g.geomout) As lat,
(addy).address As stno, (addy).streetname As street,
(addy).streettypeabbrev As styp, (addy).location As city, (addy).stateabbrev As st, (addy).zip
FROM geocode('350 Fifth Ave, New York, NY 10118') As g;
 rating | lon | lat | stno | street | styp | city | st | zip
--------+-----+-----+------+--------+------+------+----+-----
(0 rows)

So that's not great.

I noticed when loading things in I got a ton of errors like this:

Unable to open tl_2013_36123_addr.shp or tl_2013_36123_addr.SHP.
tl_2013_36123_addr.dbf: shape (.shp) or index files (.shx) can not be opened, will just import attribute data.

Is that the problem? How do I fix that?


texts: How to calculate information from overlapping cars paths on map

texts:
I have multiple geometry points that identify a car path in a time interval. Some cars can make a path that overlap another car path.

I want to calculate the car speed to have different color linestring on a map and I want to calculate the average speed for the overlapping path.

Does anyone have any idea how to select the overlapping paths and calculate the average speed?

In my Postgres database table I have these columns:
car_id , timestamp(of rilevations), latitude, longitude, geometry.

This is an example of the situation:

Green stands for 40mph.

Red stands for 80mph.

How can I calculate an average speed for the overlapping path?

The goal is to create an json file that represents the traffic state.


texts: Deploying ArcGIS Pro geoprocessing service in production

texts:
I have very little knowledge about the ArcGIS platform but we have built a geoprocessing service that calculates proximity using Python. In development, we deploy the service to server using ArcGIS Pro (by connecting to the server and deploy from the job history).
But I wonder how can I deploy to production?
I mean usually when you deploy to production, you will packaged it up in zip, jar, tar, dockers, etc. Then sometimes in old fashion way, you give this file to the site engineer to deploy.
But how about ArcGIS or more specifically ArcGIS Pro?
How can I deploy this service to production?
What kind of file I can give the site engineer?
Do we need to send developer who familiar with ArcGIS Pro to the site and bring along the python script that runs in ArcGIS Pro in a laptop (or in the server), and then deploy that to the ArcGIS server?
In ArcGIS server manager, I saw there's a button for publish service which can take a .sd file, but how can I convert my python script to .sd?


texts:
I am looking at doing a spatial query where I can find when the center of one feature falls within another. For example, when the center (centroid) of a building falls within a grid tile. I know I can do this if I create a centroid of the building polygon, but I would like to do this on the polygon layer. Is there a way to do this with ST_Within, or another function?

Thanks...


texts: non-finite value supplied by optim when I use estiStAni or fit.StVariogram

texts:
I have computed an empirical spatio-temporal variogram and I want to fit it into a model I built using vgmST. However, whenever I call fit.StVariogram I get the error "non-finite value supplied by optim".

I am still a beginner in geostatistics and I don't know if the problem is in the parameters I choose for the model I want to fit, or in choosing the stAni parameter in fit.StVariogram.
I tried calling estiStAni and I still get the same error.

Is there an easier way to automatically estimate the parameters for vgmST? (by parameters I mean the nugget, psill, range and model for both the space and time variograms), or is this error related to something else?


texts: How to import shapefile into SQL Server 2008 R2 using only SQL scripts?

texts:
I am trying to import some shapefiles into SQL Server 2008 R2. I found this useful tool called Shape2SQL however I am not sure how I can import shapefile directly using only SQL scripts as this is important for build automation and deployment. So far I haven't found anything useful just yet.

Would anyone be able to point some directions?


texts: Converting files with ogr2ogr and dropping ALL attribute fields

texts:
I want to add a gpx file to a PostgreSQL/PostGIS table using ogr2ogr, but my table does not have any of the fields that are attributes of the GPX file. There are previous questions here dealing with selecting some attributes to copy using -select or -sql rather than copying all of them, but I want to copy no attributes, as any attribute fields that try to go into the table will result in an error. I know this must be possible, given the option in QGIS to "Skip attribute creation" on file conversion. How do you do it in the command line? Something like:

ogr2ogr -update -append -f PostgreSQL PG:dbname=mydb -nln mytable filename.gpx layer -select [none? 0? null? none of these work]


texts: One of my ID numbers for my polygon is not being shown when I label my layer. How do i fix it?

texts:
I have 20 polygons and when I tried to label the layer with the ID numbers only 19 appeared. For some reason my first polygons ID number is not being shown and I don't know why.

Am I missing something vitally important?


texts:
Is there a way to strip characters from a string using python?
The string is varied in length...

Examples of strings:
Interstate 10
Auxiliary road 24
State highway 14

I am looking for just the numbers.
10
24
14

I am attempting to automate labeling roads by populating symbols with the numbers stripped from the strings.

The info I am working with is in an attribute table
I am using ArcMap 10.2.

I am new to python and very interested to learn if this sort of functionality is available.


texts: How can I increase performance on inserting features with polygon geometry?

texts:
I have a problem when inserting features with polygon feature class in FGDB. This is my code:

    ... initialising some fields (does not affect performance)
    m_featureBuffer.Shape = CreateShape(feature);
    m_insertCursor.InsertFeature(m_featureBuffer);

CreateShape:

    protected override ESRI.ArcGIS.Geometry.IGeometry CreateShape(IAbstractFeature feature)
    {
        IArea area = (IArea)feature;

        ESRI.ArcGIS.Geometry.IPointCollection polygon = new ESRI.ArcGIS.Geometry.PolygonClass();
        foreach (Point2D point2d in area.Polygon.OuterBoundaryIs.Coordinates)
        {
            ESRI.ArcGIS.Geometry.IPoint point = new ESRI.ArcGIS.Geometry.PointClass();
            point.PutCoords(point2d.X, point2d.Y);
            polygon.AddPoint(point);
        }
        return (ESRI.ArcGIS.Geometry.IGeometry)polygon;
    }

In this case inserting ~1500 polygons takes about ~110 seconds. If I comment this line

m_featureBuffer.Shape = CreateShape(feature);

the insert takes ~10 seconds! If I create empty polygons and pass them in the Shape property of the feature buffer, it takes ~10 seconds too.

I tried to use the IFeatureClassLoad according to the ESRI example, but it doesn't help with the performance.


texts: Merge multiple ASCII files into single GeoTIFF file with multiple bands

texts: Is there a way to combine multiple ASCII (text) files into a single TIFF file with multiple bands (Band1, Band2 ...) using GDAL or Python?


texts: ST_CreateFishnet inside given rectangles: Getting two geometry's instead of one

texts:
My goal is to get a fishnet inside the rectangles of 'ov_boxes' (image below) with the given extends from the attribute-table ('width' and 'height' in relation of 'scale').

I used ST_CreateFishnet and modified it:
But the function produces two geometry_columns: One have got the right spatial placement but isn't a fishnet (left image). The other is a fishnet but with the wrong (no) spatial placement (right image).

The query I used:

CREATE OR REPLACE FUNCTION ST_CreateFishnet(
        nrow integer, ncol integer,
        xsize float8, ysize float8,
        x0 float8 DEFAULT 0, y0 float8 DEFAULT 0,
        OUT "row" integer, OUT "col" integer,
        OUT geom geometry)
    RETURNS SETOF record AS
$$
SELECT i + 1 AS row, j + 1 AS col, ST_Translate(cell, j * $3 + $5, i * $4 + $6) AS geom
FROM generate_series(0, $1 - 1) AS i,
     generate_series(0, $2 - 1) AS j,
(
SELECT ('POLYGON((0 0, 0 '||$4||', '||$3||' '||$4||', '||$3||' 0,0 0))')::geometry AS cell
) AS foo;
$$ LANGUAGE sql IMMUTABLE STRICT;

drop table grid;
create table grid as
SELECT name, scale, new_geom, (ST_CreateFishnet(6,6,((width * (scale * 0.1)))/6,((height * (scale * 0.1)))/6)).*
FROM ov_boxes;
ALTER TABLE grid ADD gid serial PRIMARY KEY;

It seems to me, I produced two geometry's where I should have one: The grid inside the rectangles. I'm not very familiar with SQL, so maybe there is a miss-thinking how I create the "ouputs". But also if there is a totally different way to get to my fishnets - every idea is very welcome!


texts: How do I programmatically export a single Data Driven Page to JPEG?

texts:
I am trying to write a python script to export a individual pages to a jpeg. I am getting the error "PageLayoutObject: Error in SetDataDrivenPageID". Can anyone see where I am going wrong?

import arcpy
import os.path
arcpy.env.overwriteOutput = True
mxd = arcpy.mapping.MapDocument("CURRENT")

ref_num = arcpy.GetParameterAsText(0)
ref_num_str = str(ref_num)
ref_num_int = int(ref_num)

jobs = r"\\MSUKSERVER\BusinessMan Docs\Jobs"
job_fol = os.path.join(jobs, ref_num_str)
output_jpeg = os.path.join(job_fol, ref_num_str + "_BMS_plan.jpg")

mxd.dataDrivenPages.currentPageID = ref_num_int
arcpy.mapping.ExportToJPEG(mxd, output_jpeg, "PAGE_LAYOUT", 4679, 3308, 400)
arcpy.AddMessage("exported " + ref_num_str + "_BMS_plan")
del mxd


texts: Reproject big number of small geotiffs using gdal

texts:
I have a large number of geotiffs (a few thousand), that are in a few different projections.  IE    EPSG:32628, EPSG:32629, EPSG:32630.  I would like to reproject the tiffs to 3857 so they are all in the same projection.  I am also serving them out as 3857 so I would save some time not reprojecting on the fly.

If I reproject each file on its own gdal is very fast.  However the tiffs do not mosaic together as the reproject leaves lines at the edges.

From the image above you can clearly see where the original tiffs were stitched together.

I have thought about using gdalwarp to create one giant image, but the processing might take days and I will likely end up with a tiff too big to handle.

Is there a way to reproject the individual tiffs but take into account surrounding tiffs when doing the reprojection so I don't end up with blank lines at the edges?


texts: How to Identify and Delete duplicate vertices?

texts: I am trying to figure out how to delete duplicate vertices.  I have successfully determined that there are duplicate vertices using PLTS Data Reviewer, but I am confounded as to how I identify and remove individual vertices.  I have tried repair geometry to no avail, and I was wondering if anyone had some advice.  I am working in a file geodatabase in ArcMap 10.0 with Data Reviewer extension.


texts:
I'm attempting to read an hourly precipitation KML file from https://www.ncei.noaa.gov/metadata/geoportal/rest/metadata/item/gov.noaa.ncdc:C00313/html
in Google Collab but I am coming up short on finding libraries accessible in Collab for reading KML files. Does anyone have experience or insight into whether or not this is possible and how to do it if it is?


texts: With arcpy, how do I detect relationship classes in a file GDB (or can't I)?

texts: I would like to be able to detect relationship classes in a file GDB.  Because users of my script might have only an ArcView level license, they would not be able to manipulate the schema of a feature class (specifically to add a field) that is in a workspace that has relationship classes.  How do I detect the presence of relationship classes, so that I can document them, avoid them programatically, and allow the script to continue?


texts: How to cut polygons without lines in ArcGIS 10.4.1?

texts:
I have a road network composed of polygons and I would like to cut each polygon in way that it keeps the polygons from which the roads are made but the interior polygons are cut.
I would like to this somehow in an automate way, because doing it manually with the cut polygons tools from the edit toolbar it would take me a lot of time and I could miss some of the smaller polygons.
I have this:

and I want to cut it like this:

Lets say I have hundreds of this, all different sizes.


texts: Errno 13 "Permission Denied" when running tools and plugins on QGIS

texts:
I am having an issue with multiple tools when using QGIS where I am constantly getting PermissionError: [Errno 13] Permission denied type errors.
A lot of it seems to be associated with tools that use the GRASS or SAGA libraries, but it is also happening with some of the core QGIS Python scripts.
The example below is the error output from when I tried to run a plugin tool that calculates flood depth:
Traceback (most recent call last):
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Algorithm.py", line 409, in processAlgorithm
Grass7Utils.startGrassSession()
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 462, in startGrassSession
Grass7Utils.createTempMapset()
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 301, in createTempMapset
folder = Grass7Utils.grassMapsetFolder()
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 277, in grassMapsetFolder
folder = os.path.join(Grass7Utils.grassDataFolder(), 'temp_location')
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/algs/grass7/Grass7Utils.py", line 288, in grassDataFolder
mkdir(tempfolder)
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 103, in mkdir
mkdir(head)
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 103, in mkdir
mkdir(head)
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 103, in mkdir
mkdir(head)
[Previous line repeated 1 more time]
File "/Applications/QGIS3.10.app/Contents/MacOS/../Resources/python/plugins/processing/tools/system.py", line 105, in mkdir
os.mkdir(newdir)
PermissionError: [Errno 13] Permission denied: '/var/folders/pn'

There were errors executing the algorithm.
There were errors executing the algorithm.
Execution failed after 2.85 seconds

How do I get past this "Permission Denied" error? None of the stuff I have found so far online has provided a clear solution to the problem.


texts: Where do I find the variable name to use in a CQL statement

texts:
For example I have been following tutorials and see examples such as
    DISJOINT(the_geom, POLYGON((-90 40, -90 45, -60 45, -60 40, -90 40)))

I have accessed some map data from sources such as data.gov.uk and loaded them into an instance of geoserver where it gives the option to run a CQL statement against the map, and I have run some very basic statements to confirm it actually works.

Some data sets are point data others are polygon data, how do I find out what to use in place of "the_geom" in the example above for any given map layer I have?

PS: link to tutorial http://docs.geoserver.org/latest/en/user/tutorials/cql/cql_tutorial.html


texts: GRASS r.in.lidar (las to raster) inside Python loop?

texts:
I have a loop in Python to list out LAS files in a folder, and run two commands in Grass GIS on these files to write them out as rasters. The loop runs but at the end I get an arg must be string not list error. Please have a look at the code and suggest a work around.

import glob
import os

os.chdir(r'C:\Program Files\QGIS 3.0\bin')
grassbat = r'C:\Program Files\QGIS 3.0\bin\grass74.bat' #not needed as the os dir is changed
input_las = r'E:\Rohan\Heimdall\data\LAS'
output_directory = r'E:\Rohan\Heimdall\data\Raster\GRASS'
grass_location = r'E:\Rohan\Heimdall\grassdata\lidar\PERMANENT' # r'  reads a path as a raw string
print(grass_location)
path = "E:/Rohan/Heimdall/data/LAS/*.las"
files = glob.glob(path)
print(files)  # okay. But lists files with odd slash pattern e.g'E:/Rohan/Heimdall/data/LAS\\PNOA_2016_CAT_huso31_346-4562.las'

def las_to_GTIFF_lidardata(input_las, grass_location, output_directory):
    las_file_list = glob.glob(input_las + '\\*.las')  # lists all files in folder that (*.las)
    print(las_file_list)
    for files in las_file_list:
        tif_file_name = os.path.basename(files)[:-4] + '.tif'
        out_tif = os.path.join(output_directory, tif_file_name)
        las2raster_commands=[
            [r'grass74 -text -f E:\Rohan\Heimdall\grassdata\lidar\PERMANENT --exec r.in.lidar -e -n -o --overwrite --verbose input= files output=lidar_in resolution=2'], # removed + files here as only str is taken not list
            [r'grass74 -text -f E:\Rohan\Heimdall\grassdata\lidar\PERMANENT --exec r.out.gdal --overwrite --verbose -m input=lidar_in@PERMANENT output=out_tif + format=GTiff type=Float32 nodata=-9999'],#brackets needed?
        ]
    for command in las2raster_commands:
        print(command)
        command_final = ['grassbat', '-text', '-f', grass_location, '--exec'] + command
        print(command_final)
        os.system(command_final)
    return True

las_to_GTIFF_lidardata(input_las=input_las, grass_location=grass_location, output_directory=output_directory)

I also have tried using pyQGIS but without success:

Batch scripting r.in.lidar in pyQGIS for LiDAR DEM output?.
PyQGIS code for r.in.lidar in pyQGIS for LiDAR DEM output?.
Batch processing LiDAR files (.las to raster) with PyQGIS?.


texts: Is there a way to compute an Origin-Destination Matrix with multiple start point attributes?

texts: [QGIS 3.16] I am looking to compute an origin-destination matrix. However, from the starting point layer, I need to have multiple attributes within the matrix, opposed to the 1 linking attribute. Is this type of computation possible? I am not familiar with the log-side of coding the matrix.


texts: Arcgis Gis :Change Color on clicking Graphics

texts:
I have a number of point graphics on a map with different colors depecting some meaning.My query is on clicking any graphic(s) I want to highlight it but if I click again it should revert to its orginal color.
Please provide some inputs.
Env:
Arcgis javascript 3.14


texts: How to Display Pie Chart when clicking on point feature in map?

texts:
I want to display pie charts from MS SQL Server database data.

When I click on a point in the map, I wish to see a pie chart over the point.

Can you please tell me how to implement such pie charts?

Do you know any map site that using charts?


texts:
I have donut buffers generated to protect anonymity of household locations. They need to be clipped so that they are only inside the census tract that the original point was in. The code takes about 8 hours on a very fast computer.
Any optimizations to the following code would be greatly appreciated:

''' Clipbuff.py
'''
# Import arcpy module
import arcpy, sys, os, datetime

print datetime.datetime.now()
# Set up workspace and scratch workspace
workspace = "C:/Users/percy/projects/CEE/PerturbationProcessDataset/jitters.gdb/UTM/"
arcpy.env.workspace = workspace
scratch_workspace = "C:/users/percy/projects/CEE/PerturbationProcessDataset/scratch/"

# Allow overwriting
arcpy.env.overwriteOutput = True

# Local variables:
points = "OHAS_household"
id_Field = "OBJECTID"
polygons = "Geography_Tract"
buffers = "buff10"

# Open a search cursor on the buffers layer
rows = arcpy.SearchCursor(buffers)
print rows
# Create a for-loop to iterate through buffer IDs (via searchcursor),
# make feature layers of both the buffer and census tract layers,
# select containing census tract,
# clip by selected tract
x = 1
for row in rows:
    myVal = row.getValue("SAMPN")
    mySelection = '"SAMPN" IN ('+ str(row.getValue("SAMPN")) +")"
    # Abstract the current row value
    print datetime.datetime.now()
    # Make a feature layer with the buffers layer, selecting the OBJECTID of the     current row/record
    arcpy.MakeFeatureLayer_management(buffers, "buffers_layer", mySelection)

    # Make a feature layer with the corresponding points layer
    arcpy.MakeFeatureLayer_management(points, "points_layer", mySelection)

    # Make a feature layer for the census tracts layer
    arcpy.MakeFeatureLayer_management(polygons, "polygons_layer")

    # Select by location on the census tracts layer
    arcpy.SelectLayerByLocation_management("polygons_layer", "CONTAINS", "points_layer")

    # Clip the queried buffer by its spatially coincident census tract
    arcpy.Clip_analysis("buffers_layer", "polygons_layer", scratch_workspace + "buffclip" + str(x))
    # Add 1 to the file suffix
    x += 1
    # Delete the feature layers
    arcpy.Delete_management("polygons_layer")
    arcpy.Delete_management("buffers_layer")
    arcpy.Delete_management("points_layer")

# Set the scratch workspace as the current workspace to list the clipped buffer featureclasses
arcpy.env.workspace = scratch_workspace

# Merge all clipped buffers into one layer
mylist = arcpy.ListFeatureClasses("buffclip*")
mytarget = mylist[0]

del mylist[0]

arcpy.Append_management(mylist, mytarget, "NO_TEST")
print datetime.datetime.now()

New version of code:

# Import arcpy module
import arcpy, sys, os, datetime

print datetime.datetime.now()
# Set up workspace and scratch workspace
workspace = "C:/projects/CEE/PerturbationProcessDataset/jitters.gdb/UTM/"
arcpy.env.workspace = workspace
scratch_workspace = "C:/projects/CEE/PerturbationProcessDataset/scratch1.gdb/"

# Allow overwriting
arcpy.env.overwriteOutput = True

# Local variables:
points = "OHAS_household_subset"
id_Field = "OBJECTID"
polygons = "Geography_Tract"
buffers = "buff1"

arcpy.MakeFeatureLayer_management("Geography_Tract", "polygons" )
arcpy.MakeFeatureLayer_management("buff1", "buffers")
arcpy.MakeFeatureLayer_management(points, "points")

# Open a search cursor on the buffers layer
rows = arcpy.SearchCursor(points)
x = 1
#for row in rows:
row = rows.next()
myVal = row.getValue("SAMPN")
mySelection = '"SAMPN" IN ('+ str(row.getValue("SAMPN")) +")"
print mySelection
#print row[0]
print datetime.datetime.now()
arcpy.SelectLayerByAttribute_management("points", "NEW_SELECTION",mySelection)
ptcount = arcpy.GetCount_management("points")
arcpy.SelectLayerByLocation_management("polygons", "CONTAINS", "points")
arcpy.SelectLayerByAttribute_management("buffers", "NEW_SELECTION",mySelection)
polycount = arcpy.GetCount_management("polygons")
buffcount = arcpy.GetCount_management("buffers")
print ptcount, polycount, buffcount
arcpy.Clip_analysis(polygons, buffers, scratch_workspace + "clippy" + str(x))


texts:
This is very intriging, I don't remember having seen this before.
I am working with the Vector -> Geoprocessing Tools -> Buffer tool in QGIS (2.8.6-Wien)
If operating in British National Grid projection, when buffering a road (the road is also EPSG:27700) using 5000m, I get 5km on each side of the road when measuring the buffer distance with the Measure line tool.
If operating in World Mercator, when buffering the same road but re-projected first (EPSG:3395), using the same 5000m, I get 2.5km on each side of the road, 5km in total. It shouldn't be 10km in total again?
Why the buffer tool is using the distance buffer differently, depending on the projection. and which one is incorrect?


texts: Is there a good example of using GDAL_WMS to add an authenticated TMS layer to QGIS?

texts:
I've been trying to add a TMS layer in QGIS using a GDAL_WMS xml configuration file.

I added UserPwd to the xml and this seems to provide the authentication required by GeoServer but I'm struggling to get it to display correctly and in the correct place.

The map mostly loads as black but occasionally I can find the odd area that displays bits that look a bit like a map however I've yet to get it to load even close to correctly in QGIS.

Does anyone have a specific example of how they successfully added a GeoServer TMS layer using this method?

Example XML

<GDAL_WMS>
<Service name="TMS">  <ServerUrl>http://wms.domain.com/gwc/service/tms/1.0.0/workspace%3Alayername/${z}/${x}/${y}.png</ServerUrl>
</Service>
<DataWindow>
    <UpperLeftX>-20037508.34</UpperLeftX>
    <UpperLeftY>20037508.34</UpperLeftY>
    <LowerRightX>20037508.34</LowerRightX>
    <LowerRightY>-20037508.34</LowerRightY>
    <TileLevel>22</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>top</YOrigin>
</DataWindow>
<Projection>EPSG:3857</Projection>
<BlockSizeX>256</BlockSizeX>
<BlockSizeY>256</BlockSizeY>
<BandsCount>3</BandsCount>
<Cache />
<OfflineMode>false</OfflineMode>
<UserPwd>user:password</UserPwd>
<ZeroBlockHttpCodes>204,303,400,404,500,501</ZeroBlockHttpCodes>
<ZeroBlockOnServerException>true</ZeroBlockOnServerException>


texts: Seeking qml or sld file for QGIS + OpenStreetMap data?

texts:
I'm looking for a qml file (or equivalent e.g. sld) to allow me to use QGIS to take OpenStreetMap data and produce paper maps. I'm re-asking this question (it was asked once or twice over a couple of years by other people) because I'm really surprised not to easily find such a thing already (I looked hard).

I have found Anita Graser's styles - which are good, but only style basic data like roads (and do so simply). I've successfully followed the directions for getting OSM data to work with. I've carried out some basic style editing successfully. I realise I could create my own styles for this purpose - but this seems like a very big job. I want a visually pleasing map taking into account most of the data available on OSM (i.e. not just streets, but tracks, paths, rivers, woods, lakes, buildings and so on). I'd have thought that there would be other people out there already doing the same... It seems odd to have to re-invent the wheel when both OSM and QGIS are open source / open data projects. I feel like maybe I've just missed the obvious repository of OSM qml files which resides out there on the internet somewhere if you know where to look.

Additional notes in response to query below...

Ideally I'm looking for a Mapnik style rendering of the OSM data, but any equivalent properly developed and complete style would be good.

At the moment the objective is to use this data to produce paper mapping, but simple use of the mapping as base data on which other geographic data can be displayed also requires a decent rendering. At work I use Ordnance Survey data completely styled to produce a full and very detailed UK map (the point being that I didn't have to design this styling, it came with the data). I realise that other tools exist for paper mapping, but the same issues tend to arise at one stage or another... I'm more or less successful in making the tool work for me in really basic terms, but to get to something visually useful I end up back with me either needing to re-invent a complete rendering style from scratch (beyond my skill), or encountering such a complex set of requirements for database setups or whatever that the process crashes and burns before I get anywhere useful. With a few more years of GIS experience I'll maybe be a database, CartoCSS (or whatever else) wizz... but until then I'm prevented from using the amazing resource that OSM provides by this one barrier.

Update (July 2015): Please note that I'm still looking for information on this, and given that I've just been informed that the question has had 2500 views in the last year clearly others are too.


texts: Load coordinates from database to Google Maps

texts:
How can I load coordinates (lng/lat) from a database (MySQL) into a map service like Google Maps or Leaflet?

Currently I only have a JavaScript which loads a GeoJSON file with coordinates.

I want to inject them from my database to the map. I'm not very skilled in web development.

Why can't I just do SQL queries in JavaScript and print the results to the map like in Java where I just loop the result set and operate on my data.

Or do I have to use PHP or some other server side scripting language? The issue is, I have to use Java to process the coordinates.


texts: ArcGIS Online HTML popup: Why does the code for the table background colors keep getting removed?

texts:
I have the following code as my popup

  {counties_NAME}

  Income spent on Housing &amp; Transportation...

{ht_index_data_v3_csv_ht_ami}%

It works fine for a while, but then all the background colors will disappear. The code then looks like this without me changing anything.

  {counties_NAME}

  Income spent on Housing &amp; Transportation...

{ht_index_data_v3_csv_ht_ami}%

ArcGIS online is removing the style portions of my code and I'm not sure why. Online, it says it's supported.


texts:
I am having a problem georeferencing. I have a photograph (.tif) that will not save georeferencing. When the .tif raster appears on the table of contents, there is a ":0" added to the end of the naming convention. I then fit to display and start adding control points. I can create a text from link table, and rectify file no problem. The problem arises when I go to update georeferencing. The image will disappear along with the control points. When I search for the raster, I find it @ 0,0. Any ideas?

I found a similar post from a few years back. Same problem but his solution didn't help my issue...


texts: Classify values in attribute table that is joined to a raster in ArcGIS 10?

texts:
I have a an attribute table that does not have a raster value and is joined to another table that does. How do I reclassify the values in the attribute table?

For example, in column B, I'd like all values <5 to reflect the value 1. How do I do this classification?


texts: Updating multiple feature classes from one table

texts:
I have a large citywide XL table with address and meter information that I broke out spatially by neighborhood association into 29 different feature classes to make usage maps. It was a fairly long process (geocoding, several spatial joins, etc) to get things all broken out.

Now I have an updated version of the original table with updated meter reads, and I'm looking for a way to automate populating the new meter reads into the individual feature classes. This may end up being a monthly job. I have a field I can join by, but I'm stuck as to how I can do this across 29 feature classes without repeating the process 29 times.

Is there a way I can use ModelBuilder or (fairly basic) python to run this process across multiple separate feature classes in one fell swoop?


texts: Filtering features in other layers of QGIS Atlas?

texts:
I'm trying to use the QGIS Atlas tools to build student 'scattermaps', ie. display the school and the students that attend that school, for 215 schools.

The schools layer is the coverage layer, and is being driven by the SCHNUM (school number) value. This is working fine, that is it is driving the pages correctly.

How, though, can I display only features from other layers, in this case the student point layer, that match the corresponding SCHNUM value for the current feature in the coverage layer?

I would also like to switch the coverage layer (schools) to 'hidden' and use the same technique to only display the current school from a duplicate schools layer.

There is a presentation available here that alludes to the use of a layer rendering rule setting using $atlasfeatureid, but I can't figure out how that works.

I also found this bug tracker page that shows this working, but again, I can't get it to work on my data:


texts:
Can I project a NetCDF file with the help of GDAL or QGis?
The file info is:

gdalinfo c:\tomcat-6.0.36\webapps\ncWMS\layers\RAIN_MM_1hour_ac.nc
Driver: netCDF/Network Common Data Format
Files: c:\tomcat-6.0.36\webapps\ncWMS\layers\RAIN_MM_1hour_ac.nc
Size is 772, 668
Coordinate System is `'
Origin = (19.980464188539599,49.004338539807151)
Pixel Size = (0.012931622920800,-0.008983060104580)
Metadata:
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.5
  NC_GLOBAL#GDAL=GDAL 1.11.2, released 2015/02/10
  NETCDF_DIM_EXTRA={Time}
  NETCDF_DIM_Time_DEF={1,6}
  NETCDF_DIM_Time_VALUES=0
  RAIN_MM_1hour_ac#_FillValue=0
  RAIN_MM_1hour_ac#long_name=GDAL Band Number 1
  RAIN_MM_1hour_ac#units=mm
  Time#units=hours since 2014-10-01 00:00:00
Corner Coordinates:
Upper Left  (  19.9804642,  49.0043385)
Lower Left  (  19.9804642,  43.0036544)
Upper Right (  29.9636771,  49.0043385)
Lower Right (  29.9636771,  43.0036544)
Center      (  24.9720706,  46.0039965)
Band 1 Block=772x1 Type=Float64, ColorInterp=Undefined
  NoData Value=0
  Metadata:
    _FillValue=0
    long_name=GDAL Band Number 1
    NETCDF_DIM_Time=0
    NETCDF_VARNAME=RAIN_MM_1hour_ac
    units=mm

Using gdalwarp

gdalwarp -t_srs '+proj=lonlat +datum=WGS84 +no_defs'
 c:/tomcat-6.0.36/webapps/ncWMS/layers/RAIN_MM_1hour_ac.nc
 c:/tomcat-6.0.36/webapps/ncWMS/layers/RAIN_MM_1hour_ac_proj.nc

results in:

ERROR 1: Failed to lookup UOM CODE 0
Warning: The target file has a 'nc' extension,
which is normally used by the GMT, netCDF drivers,
but the requested output driver is GTiff. Is it really what you want ?

I tried to save the file in QGis with the specified projection (EPSG:4326), but the output file is limited to a GeoTiff format.

I want the output data to also be in NetCDF format, in order to serve it to an OpenLayers map through ncWMS.

PS: I am not familiar with R or Python yet, so I would like an approach using the above mentioned software.


texts: PyQgis on stand-alone Python software says ImportError: No module named qgis.core?

texts:
I'm quite new to Python, but I learned the basics with ArcPy in the Python software IDLE. It was quite easy to getting started: I start my script with "import arcpy" and i could start coding with arcpy functions.

But I'd like to work with QGis and PyQgis now. And I cannot make it work. When I start my script with "from qgis.core import *", i get an error message saying:

Traceback (most recent call last):
  File "<pyshell#1>", line 1, in <module>
    from qgis.core import *
ImportError: No module named qgis.core

I guess that I have to make a liaison between Python and the Qgis module somehow, it is probably very basic but I can't find the answer.

How do I set-up my python so it can use the qgis functions?