texts:
mapproxy.config - WARNING - unknow 'html' in sources.gasps_osm_source has thrown by terminal
when i run : mapproxy-seed -c 16 --continue -f mapproxy.yaml -s seed.yaml , and no one tile has been cashed... o tiles.

settings:

services:
  demo:
  tms:
    use_grid_names: true
    origin: 'nw'
  kml:
      use_grid_names: true
  wmts:
    # use restful access to WMTS
    restful: true
    # this is the default template for MapProxy
    restful_template: '/{Layer}/{TileMatrixSet}/{TileMatrix}/{TileCol}/{TileRow}.{Format}'
    # and also allow KVP requests
    kvp: true
  wms:
    md:
      title: MapProxy WMS Proxy
      abstract: This is a minimal MapProxy example.

layers:
  - name: osm_background
    title: osm_background
    sources: [gasps_osm_cache]

caches:
  gasps_osm_cache:
    grids: [osm_grid]
    sources: [gasps_osm_source]
    meta_size: [6, 6]
    meta_buffer: 100

sources:
  gasps_osm_source:
    type: wms
    html:
      client_timeout: 500
    req:
      url: 'http://xxx.local:8080/geoserver/test_osm/wms'
      layers: osm:osm_layerset
    coverage:

      srs: 'EPSG:3857'
grids:
    osm_grid:
        base: GLOBAL_MERCATOR
        srs: 'EPSG:3857'
        origin: nw

globals:

and

seeds:
  myseed1:
    caches: [gasps_osm_cache]
    # grids: []
    #coverages: [changed_area]
    levels:
      to: 12
    refresh_before:
      time: 2019-9-25T12:35:00


texts: Sum Grouped by Query of Values from Field of Variable Length

texts:
I am trying to find the sum of watershed areas for about 30 small streams in BC.  Each stream has a code that is the same as the code for its watershed.  The code for each stream is derived by adding several digits to the code of it's parent stream.  For example; stream [c] 100-123-456-789-000-000-000 is a tributary of stream [b] 100-123-456-000-000-000-000, which is a tributary of stream [a] 100-123-000-000-000-000-000.

So, if streams [b] and [c] are the only tributaries of stream [a], then the total area of watershed [a] is the sum of the area of all three streams.

I think that the best way to perform this query would be using SQL in the Field Calculator (QGIS).  I am trying to write the sql query myself, but I am having trouble grouping the streams because of the varying length of the code using:

I think that I would need to use something like replace to deal with the varying length of code. So I create a virtual field (code.trunc') with a truncated version ofcode' using : 'Replace("code", ’-000’,’’).

Then I am trying to sum the watershed area of the tributaries of each stream using:

sum ("Sum.Val", "code.trunc" LIKE concatenate("code.trunc",'%')

This code crashes QGIS.  I am guessing that it is because the LIKE function is not very efficient, and this is a large data set.  Note, I am using a subset (n=6144 features) of the data set which I will link to below).

Is there a better way to do this?

This table shows what the results of the query might look like:

The data that I am using is available here.  When I figure out the query that I need to use, I plan to perform it only on the 30 streams which will be selected prior to the query.


texts: Error with the Contour (Spatial Analyst) tool in ArcGIS

texts:
I am having an issue with the Contour (Spatial Analyst) tool in ArcGIS. When I manually run the tool in ArcToolbox, it will sometimes say that the tool has failed, giving the following error message:

ExecuteError: ERROR 010157: Unable to open feature class F:\2014\data\083c15_contour.shp
Failed to execute (Contour).

The strange part is that when I go to the output location of the tool the shapefile is there, and I can open it in ArcMap without an issue.

When I use Python to run the tool, I receive the same error message, and I can't figure out why. This is the loop I have that is meant to convert a list of DEM files into contour shapefiles:

    for dem_file in dem_files:
        file_split = dem_file.split("_")
        file_name = job_folder + "\\" + file_split[0] + "_contour.shp"
        arcpy.sa.Contour(dem_file, file_name, contour_interval)

When I run the script, I can go to the output location and bring the shapefile into ArcMap, even though there was an error. So I don't think this is an issue with my code because it also happens when I run the tool manually. Does anyone know how I can work around this error?


texts: Exporting current map extent or bookmark to polygon feature using existing tool?

texts:
Is there a out-of-the-box tool in ArcGIS 10 (Editor) that would allow me to easily create a polygon feature from the current map extent or to convert a bookmark to a polygon feature.

I am not interested in a programmatic or a 3rd party solution.


texts:
I am quite new to QGIS and the Time Manager Plugin.
I am currently searching for a way how I can show satellite data images (sentinel 2 for example) with the time manager plugin.
The idea is to have an image of a certain area every day/month/year and with the time manager I want to show the changes in terms of clouds, growing buildings, or movement of ships for example.

Is that possible?
Currently, I only saw examples of the time manager plugin where shapefiles and vector data got used.


texts: Bring marker on same layer to front in CartoDB

texts:
I'm using a category map showing markers based on the value of one column. They can be one of two values. One is for customer location and the other for store location.

I'd like to have it so that the store locations are ALWAYS shown on top of the customer locations. Is there a way to force this with CartoCSS. Looking for something like z-index but this doesn't seem to work.


texts:
I came by an "infiltration map" today.

But I cannot find any definition of this on Wikipedia or the web.

What is it?


texts: Iterate through features of a feature class and buffer them?

texts:
I'm currently trying to iterate through the features of a feature class (using an update cursor, because the goal is to ultimately update the field) and buffer each individual segment. I eventually want to use the buffer in each iteration with different tools, but I can't get the initial Select by Attribute to work. Here's my code:

fc = "Original_Clip"
fcfield = "SEGMENTID"

strt_lyr = arcpy.MakeFeatureLayer_management(fc, "Street_lyr")
with arcpy.da.UpdateCursor(fc, fcfield) as cursor:
     for row in cursor:
        var_segment = row[0]
        print "Buffering segment " + str(var_segment)
        arcpy.SelectLayerByAttribute_management(strt_lyr, "NEW_SELECTION", '{} = \'{}\''.format(var_segment))
        arcpy.Buffer_analysis(strt_lyr, str(var_segment) + "Buff", "20 Feet")

When I run this code, an error is thrown with the SelectLayerByAttributefunction, saying my where clause is invalid?

Runtime error Traceback (most recent call last): File "", line 15, in File "d:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 6688, in SelectLayerByAttribute raise e ExecuteError: ERROR 000358: Invalid expression Failed to execute (SelectLayerByAttribute).
I'm assuming I just have the SQL syntax wrong (I'm querying a file geodatabase btw), but maybe what I'm trying to do just isn't possible.


texts:
In ArcMap I have a data frame with a base map using the WGS_84_Pseudo_Mercator projection.

When I add a shapefile layer with a WGS_84_Pseudo_Mercator projection, the shapefile position appears slightly off, in respect to the base map - About 200 pixels off.

I've tried defining a new projection, and even reprojecting the shapefile layer, but it always appears in the same position.

Everything seems to be correct, so I have no idea where should I be looking for - I've even tried to load the shapefile in QGIS, and it appears in the correct position.

Any ideas?


texts: Finding distance of point along polyline using ArcGIS Desktop?

texts:
We need to find the distance of a series of point along a given polyline.

We have a polyline dataset repsenting rivers in a catchment, broken up into river reaches (FID 1 through 35 for example).

We also have a point dataset of abstraction licenses within this catchment, some of which are situated along rivers - we can tell from the attribute table which these are as they are defined as groundwater or surface water licenses.

The two are not related - i.e. one has not been derived from the other.
Therefore not all the surface water licenses lie exactly over the river polylines, despite in 'reality' being on these rivers.

We need to be able to see which points are "upstream/downstream" of each other, and so wanted to have a list of which river reach the point is on, followed by the distance downstream from the start of this reach. We can then compare the downstream distances of two point to see which is upstream and which is downstream - the desired outcome!

To solve the problem of finding points on/very close to the rivers we were told to use a (small) buffer around each river reach and then to query the buffer to find the points within it.

This is fine but we are then left with an area (buffer), not a line, so we can't measure a distance along/within this.

Can we use a buffer as described to find out the points on/near the river reach in question and then, somehow (?!) query these points for distance along the river reach (polyline segment)?


texts: Count points not in polygon in QGIS (difference between points and layers give unexpected result)

texts:
I have a layer of 130 points and another layer of driving distances as 6 polygon features
When I use the count points in polygons I get a running total of the number of points in each polygon feature (as expected) with the largest number being 100

Difference

When I get the different of the same points and the polygon layer I get a layer containing more points than I expect (expect 30 but get 50)

Is difference the correct way to select points not in layer and what could be the reason for the inaccurate result (when I visualise the layer I don't see a symbol for each of the 50 points just an handful

  Both the points and polygon layers are using the same CRS


texts:
I am using QGis 2.18.7 and would like to open a jpg from ncap.org.uk:
Here is the link to the image but one has to register in order to download the image:
https://ncap.org.uk/frame/6-1-17-30-391
There is a jpg and a kml file for download. I put both in the same folder but I can not find a way to open the jpg in QGis.


texts:
So i am using leaflet in order to create my project, which is a world map that certain layers would be add on top of it.

I am wondering whether I should change the CRS property of the map object, which is by default EPSG3857. The map I am using as tilelayer is from openstreetmap, whose maps i think are in EPSG3857. However, the layers I want to add on top of the map, are in EPSG 4326 when I view them on QGIS.

For example, this is my map object:

var map = L.map('map').setView([40, 25], 9);
L.tileLayer('http://{s}.tile.openstreetmap.se/hydda/base/{z}/{x}/{y}.png', {
//crs: L.CRS.EPSG4326,
attribution: 'Tiles courtesy of OpenStreetMap Sweden &mdash; Map data &copy; OpenStreetMap',
}).addTo( map );

and I load the layers from Geoserver, like this:

var randomlayer = L.tileLayer.wms("http://localhost:8080/geoserver/wms", {
layers: 'workspace_name:random-layer',
format: 'image/png',
transparent: true,
version: '1.1.0',
attribution: "myattribution"
});

Leaflet clearly mentions about CRS "Coordinate Reference System to use. Don't change this if you're not sure what it means." and that EPSG3857 is "The most common CRS for online maps, used by almost all free and commercial tile providers. Uses Spherical Mercator projection. Set in by default in Map's crs option."

So basically, should I change the crs property of my map object, aka add:

crs: L.CRS.EPSG4326,

on it, in order for it to be the same as the layers, or should I leave it as default which I think is the same as the map tilelayer?


texts:
I saw an example of "Render geometries to a canvas":

  var canvas = document.getElementById('canvas');
  var vectorContext = ol.render.toContext(canvas.getContext('2d'), {size: [100, 100]});

  var fill = new ol.style.Fill({color: 'blue'});
  var stroke = new ol.style.Stroke({color: 'black'});
  var style = new ol.style.Style({
    fill: fill,
    stroke: stroke,
    image: new ol.style.Circle({
      radius: 10,
      fill: fill,
      stroke: stroke
    })
  });
  vectorContext.setStyle(style);

  vectorContext.drawGeometry(new ol.geom.LineString([[10, 10], [90, 90]]));
  vectorContext.drawGeometry(new ol.geom.Polygon([[[2, 2], [98, 2], [2, 98], [2, 2]]]));
  vectorContext.drawGeometry(new ol.geom.Point([88, 88]));

But what to do with the geometry in the projection EPSG:4326?


texts: Flow Mapper wrong line length calculation

texts:
I'm finding some trouble with the length calculation in the Flow Mapper plugin. As you can see in the picture below the length of the 2 lines

C001-C021 = 4233.7986...

C001-C009 = 6069.3053...

is wrong, as you can tell even just looking at the map. To confirm this I've output the distant matrix using the tool found in the vector menu -> analysis -> distance matrix that you can see in the excel spreadsheet in the top-right of the picture.

Why do these calculations give me different results?


texts: BlenderGIS no longer creates fine grid of elevations from SRTM30

texts:
On my first effort using BlenderGIS, I could get fine detailed SRTM30 elevations of this peninsula at 400,100 triangles but lately I'm only getting roughly detailed SRTM30 elevations of the same peninsula no better than 8192 triangles. Shown are the two side by side after converting to DXF.

How can I obtain the track log requesting the DEM server to see what was requested and what was downloaded?
I'm using Blender version: 2.91.2, branch: master, commit date: 2021-01-19 16:15, hash: 5be9ef417703, type: Release
build date: 2021-01-19, 16:25:50


texts:
When I run this code, it says "An invalid SQL statement was used". Any idea what does it mean?

import arcpy
from arcpy import env

env.workspace ="C:/RG/Projects/CriticalDataNetwork/CDEC.gdb"

try:
    inTable = "Avg_Daily"
    keyField = "STATION_ID"
    inField = [["Avg_Daily.OBJECTID","OBJECTID",],["Avg_Daily.STATION_ID","STATION_ID"],   ["Avg_Daily.DLY_AVG", "DLY_AVG"]]
    whereClause = "\"STATION_ID\" = 'CAP', 'DDM'"
    outTable = "Test"
    print outTable

    arcpy.MakeQueryTable_management(inTable,  outTable, "USE_KEY_FIELDS", keyField, inField, whereClause)

except Exception, e:
    # If an error occurred, print line number and error message
    import traceback, sys
    tb = sys.exc_info()[2]
    print "Line %i" % tb.tb_lineno
    print e.message


texts: Using pgsql2shp in server code with PostGIS installed on Amazon RDS instance

texts:
I've created a web application using the Meteor web framework, GeoServer and PostGIS. I've deployed my app on Amazon Web Services with the Meteor app and GeoServer running as 2 separate tasks on 1 EC2 instance using ECS, and my PostGIS database being hosted by Amazon RDS. I've run into an issue where I've added a feature allowing users to download a Shapefile of their data by calling pgsql2shp in my server code, but this only works locally for me as I have PostGIS, GeoServer and my app all running on my machine. In the deployed version of my app the download feature is causing my app to crash as I don't have PostGIS installed on the container running my server code, it's installed on the Amazon RDS instance.

Does anyone know how I can export a Shapefile from Amazon RDS using pgsql2shp, when the utility isn't installed on the server machine that is executing the export command? Is there a web API that can be called instead?


texts: Adding page to layout after populating the page using PyQGIS

texts:
TL;DR How can I add items to the page (QgsLayoutItemPage) without using layout?
I want to add an item to the 2nd (or 3rd, ...) page in the project layout. Normally, I can add the item to the 2nd page by specifying the layout coordinates using the following script in this way:

Initialize a layout,
Add new page (2nd page) and
Add items to the layout (using the layout coordinates)

import os
from qgis.PyQt.QtGui import QFont

layout_name = "LAYOUT_NAME"
project = QgsProject.instance()

layout_manager = project.layoutManager()
layouts = layout_manager.printLayouts()

# remove the same named layout
for layout in layouts:
    if layout.name() == layout_name:
        layout_manager.removeLayout(layout)

### (1) Initialize a layout
layout = QgsPrintLayout(project)
layout.setName(layout_name)
layout.initializeDefaults()
layout_manager.addLayout(layout)
iface.openLayoutDesigner(layout)

### (2) Add new page (2nd page)
page = QgsLayoutItemPage(layout)
page.setPageSize('A4', QgsLayoutItemPage.Landscape)
layout.pageCollection().addPage(page)

### (3) Add items to the layout (using 'layout')
label = QgsLayoutItemLabel(layout)
label.setText("TEST")
label.attemptMove(QgsLayoutPoint(10, 230, 0))
label.attemptResize(QgsLayoutSize(10, 10, 0))
# ADD ITEM
layout.addLayoutItem(label)
###

For example, the coordinates would be (10, 230) to add a label to the 2nd page for a A4 landscape size 297x210 mm (see the image).
230 = 210 (height of the 1st page) +
      10  (space between pages) +
      10  (y value for the label in the 2nd page)

But I want to add items using QgsLayoutItemPage object so that I don't have to deal with calculating the layout coordinates as in this answer.
I would like to follow this way:

Initialize a layout,
Create new page, don't add to the layout,
Add items to the page and
Add the page to the layout

How can I do that?


texts: Publishing ArcGIS map services to a different server

texts: Is it possible to publish a map service to a server different to the one where ArcGIS Server was installed? If so, how can this be done?


texts: Problem to clip a raster with a shapefile

texts:
I need to clip a raster layer with a vector layer. When I try this, the  following error appears:

**ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.**

What am I doing wrong?


texts: PyQgis - Adding attributes to a POSTGIS layer

texts:
After reading @underdark's answer to this question on creating a new vector layer and @Matthias Kuhn's to this one on field calculation, I wrote a script that allows me to calcultate the mean arrival-on-scene time for emergency response services per cell of a grid, even thought some of my values are null.

Steps :

Add attributes in which to calculate mean
Sum the average times of arrival per year (I'm working on 5 years)
Get the number of years for which value is not null
Divide sum by number of not-null-values

My problem is as follows :

The script works fine on OGR layers, but not on Postgres layers. I don't get any errors, but nothing happens when I try to add fields, which is my first step. So right now I have to save my postgres layer as a .shp on my computer to do the treatments I want. Which therefore don't automatically get saved in my database...

Here is the beginning of my script (the part that creates new attributes) :

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

# supply path toqgis
QgsApplication.setPrefixPath("C:\Program Files (x86)\QGIS Wien", True)

# load provider
QgsApplication.initQgis()
qgis.utils.iface

# set active layer
clayer = qgis.utils.iface.activeLayer()
provider = clayer.dataProvider()

# start editing mode
clayer.startEditing()

# add new fields
caps = provider.capabilities()

if caps & QgsVectorDataProvider.AddAttributes:
    res = provider.addAttributes([QgsField("quotient", QVariant.Int), QgsField("moy", QVariant.Int)])

clayer.updateFields()

Is there a code line to indicate I want/have rights to alter my postgres layer ?

The most similar questions I could find - this one and this one - didn't get a solution.

I should say I am new to PyQgis so if you are so kind as to answer please explain your suggestion as clearly as possible.


texts: Is there any software (with a graphical user interface) for creating Style Layer Descriptors?

texts: Is there any software that can easily create a Styled Layer Descriptor (SLD) without using a  script?


texts: Define pen width of a Line layer in PyQGIS

texts:
I have created a line shapefile layer (memory layer) created in Python Console of QGIS. Now I want to define the thickness of the line features that are being created. I am aware of the method by which the line styles/symbology can be imported from the existing '.qml' file.

self.setLayerSymbology(vlayer, "Parent_PermanentMarks.qml")

Is there any way by which symbology (line thickness) of the lines can be defined 'on the go' via python console?


texts: Gradient buffering background layer in ArcGIS Desktop?

texts:
I have a map that was created by a contractor working for my employer before I began. I am not sure what GIS software the map was created in but it has an effect I would like to be able to reproduce in ArcGIS 10.1.

The map seems to have a gradient, buffered layer on the outside of the Australia coastline layer. This layer grades out from dark grey to white. Sitting off the coast in my layout view I have the Australia with the extent area showing.

What I can't get is the gradient layer to gradient out behind the Australia extent part. In the maps I make Australia is in a separate data frame and I can't make the background of that data frame transparent.

I can't post the maps as I'm a new user but if someone can assist me I can send you directly the maps to explain better what I mean.


texts:
How to remove all polygons in QGIS from layer A which are at least partially overlapped by any polygon from layer B?
Geoprocessing tools / Difference is not an option as it will remove only overlaping parts of the polygons from layer A and not whole polygons.


texts: ArcGIS JS Modular Feature Layers? Possible?

texts:
My ultimate goal is a map which can toggle the visibility of layers on an ArcGIS map. However the data I need to use is not available as a feature service, I have to build it programatically from feature collections. I want to make AMD modules that will return FeatureLayers to a script that is in control of the map. I also want the FeatureLayers to contain events.

Does anyone have any advice on how to accomplish this? Can I setup the events and map markers without access to the map object? Maybe pass the map object to the FeatureLayer constructor? Any advice is greatly appreciated.

EDIT: This is what the code that calls my modules looks like.

require(["dojo/parser", "dojo/_base/array", "dojo/request", "dojox/layout/TableContainer", "dijit/layout/BorderContainer", "dijit/layout/ContentPane", "dijit/layout/AccordionContainer", "esri/map", "esri/layers/FeatureLayer",
             "esri/dijit/Popup", "esri/tasks/query", "dgrid/Grid", "dojo/json", "js/cbibsModule", "dojo/domReady!"],

              function (parser, array, request, TableContainer, BorderContainer, ContentPane, AccordionContainer, map, FeatureLayer, Popup, query, Grid, JSON, cbibsModule) {
                function init(){
                    _map = new map("mapDiv", {
                      basemap: "oceans",
                      center: [-77.0357, 38.7877],
                      zoom: 7
                    });
                    var cbibs = new cbibsModule();
                    var fl = cbibs.requestCBIBS(_map);
                    console.log(_map);
                }
                dojo.ready(init);
              }
            )

cbibs.requestCBIBS(_map) should not only return a Feature Layer object, but should also pass the _map variable by reference so the Feature Layer can attach events to it. This is not working. I don't get any errors but the map does not display.

UPDATE: This is my updated code but my markers are still not drawing on the map, nor am I getting any errors. I build the Graphic objects in a module and return them in an array. I then add them to the Graphics Layer and add it to the Map (at least I think that's what's going on). Any suggestions??

          function (parser, ready, array, request, TableContainer, BorderContainer, ContentPane, AccordionContainer, map, GraphicsLayer, Popup, query, Grid, JSON, cbibsGfxModule, crwoModule) {
            var Map, gfxLayer;
            function init(){
                Map = new map("mapDiv", {
                  basemap: "oceans",
                  center: [-77.0357, 38.7877],
                  zoom: 7
                }),
                dojo.connect(Map, "onLoad", addGraphics);
            };

            function addGraphics(){
                gfxLayer = new esri.layers.GraphicsLayer();
                var gfxContainer = []; // Store array returned by cbibsGfxModule.getAllCurrentReadings()
                var cbibs = new cbibsGfxModule();
                gfxContainer = cbibs.getAllCurrentReadings();
                array.forEach (gfxContainer, function (item, i) {
                   gfxLayer.graphics.push(item);
                });
                Map.addLayer(gfxLayer);
                console.log(Map);
            };
          dojo.ready(init);
          }


texts: If I wanted to advertise for a freelancer to do a GIS job, where would be a good place? I've advertised on oDesk (which we use frequently), but it doesn't even have a GIS category. Am I allowed to do it here?


texts: Refering to the input-layer of a QGIS model in field calculator expression

texts:
The setting
In a QGIS model, I have an input layer polygon from which a grid is cretated.
On this grid-layer, a new attribute field should be created based on an expression that refers to the input layer polygon. That works well as long as I set a fixed name for the layer polygon - see the expression and the model linked here (at the bottom).
Basically, I just have to replace 'polygon' in the overlay_intersects-part of the expression in a way to refer to the model's input layer: overlay_intersects('polygon',$geometry) works, as long as the layer indeed is named 'polygon', but the model should be flexible to take the selected input layer here.
The question
How can the layer-name be defined so that it refers to the name of the input layer, regardless of how it's name is?
What I tried
In the field-calculator dialog, I select pre-defined value and try to set my expression in 'single quotes' as described here in the style of: 'expression' || @loadpolygonlayer || 'expression'. For the layer to refer, I select it from the context variables as can be seen on the screenshot.
I also tried to get the name of the layer by using this syntax and different variations with and without quotes, with masking single quotes etc., with no success.


texts: I have about 4000 stores all across the country. Based upon the long/lat distance and business level, I wish to create about 200 clusters who can manage 15-25 stores each. I wish to find out how to go about this to find these 200 clusters. Is it possible to do this? If any suggestions, please suggest how to go about it.


texts: Convert total area calculation from ogrinfo (cli) to python ogr

texts:
How can I do this command (source):

ogrinfo -q -sql "SELECT SUM(OGR_GEOM_AREA) AS TOTAL_AREA FROM way" way.geojson

with Python ogr?

from osgeo import ogr


texts: Understanding How ArcGIS Displays Non Spatial Data Created in Microsoft Access

texts:
I am working on getting American Community Survey(ACS) data into SDE table format. I have managed to import all the text files downloaded from ACS into excel and then imported to Access for further extracting.

Anyway, so all the tables have an ID field. When I try to preview the data in ArcCatalog it would not show the data in proper order shown below. The problem is when I register this tables as a geodatabase tables it is not using an ID field to create OID field but it assigns OID to the order that it was displayed in preview tab. I want to understand how it creates and assigns the OID to each row. I know I can right click on the field and sort by ascending order and it would display how I want it but it is a temporary solution. This does not solve the problem when I convert the access table to a SDE table.

Thank you in advance.


texts: Shapefile CRS changed after exporting from QGIS. Are they equivalent?

texts:
I added a shapefile as a vector layer in qGIS 1.7 so that I could join in additional attribute data (not geodata) from another DBF file.  The Shapefile that I added has this CRS (as far as I can tell using ogrinfo):

PROJCS["WGS_1984_Web_Mercator",
GEOGCS["GCS_WGS_1984_Major_Auxiliary_Sphere",
    DATUM["WGS_1984_Major_Auxiliary_Sphere",
        SPHEROID["WGS_1984_Major_Auxiliary_Sphere",6378137.0,0.0]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]],
PROJECTION["Mercator_1SP"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["latitude_of_origin",0.0],
UNIT["Meter",1.0]]

After joining in the extra attribute data in, I exported the original shapefile by right-clicking on the layer and choosing 'Save as.'  I kept the CRS as 'Original CRS' in the 'Save vector layer as ...' dialog.  However, the exported shapefile has the following CRS defined:

PROJCS["Mercator",
GEOGCS["GCS_unnamed ellipse",
    DATUM["unknown",
        SPHEROID["Unknown",6378137,0]],
    PRIMEM["Greenwich",0],
    UNIT["Degree",0.017453292519943295]],
PROJECTION["Mercator_1SP"],
PARAMETER["central_meridian",0],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["Meter",1],
PARAMETER["latitude_of_origin",0.0]]

Are the CRS' equivalent?  If so, does that mean that some of the strings in the WKS definition of the PRJ file are just identifiers?  If not, how can I make qGIS preserve the CRS of the original shapefile when exporting it?


texts: Reducing processing time of multiband satellite images using RSToolbox in R

texts:
I have created a script in R using the RSToolbox and raster packages in R to process multiband images from Landsat 7 and 8. The code works fine and it is giving me the outputs I need; however, due to the amount of data contained in these images, it takes a long time to run (~30-40 minutes per image and I have nearly 200 images to analyse).

To reduce the analysis time, I tried cropping the image to just the area I am interested in, which worked really well, but after a couple of calculations, the image becomes over 1GB in size and an error comes up saying that R is unable to process an image that size. I then noticed that when using the crop function, the image becomes about 10 times larger than the original one, even though it is cropped to about half its original spatial extent. I am not really sure why this is the case.

The image was cropped using the following code:

library(raster)
cropMB = crop(MB,RioDoce)

where MB is the stacked layer of a Landsat 7 or 8 image, and RioDoce is a point vector representing the path of a river (uploaded to R as a shapefile).

Landsat images were stacked using the following code:

library(RSToolbox)
meta = readMeta("filepath to Landsat image metadata")
MB = stackMeta(meta)

Would anyone here be able to suggest a way to decrease the time of the analysis by decreasing the size of the original image? For reference, I am doing the following with each set of images:
- Create a raster stack containing all bands
- Transform green and NIR bands from DN to TOA
- Calculate NDWI using the previously mentioned bands
- Apply two masks to original raster stack so only areas with water bodies and no clouds remain
- Transformation of DN values to at-surface-reflectance of the masked raster stack
- Extract spectral band values to a point layer representing the path of a river

How can I reduce the analysis time?


texts: TROPOMI NO2 and SO2 too small ... being approximated to zero

texts:
I am trying to download TROPOMI NO2 and SO2 data for my region as a timeseries from Google Earth Engine, however, the values are small (see attached images) and as such, when I download the CSV file, I find that all the values have been reduced to zero (0).
How can I solve this in order to download the values as they are?

Here is my code: link
var geometry = ee.Geometry.Polygon(
  [[[27.671093749999986, 5.861609095403631],
  [27.671093749999986, -12.999528581977067],
  [44.19453124999998, -12.999528581977067],
  [44.19453124999998, 5.861609095403631]]], null, false);

var data = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')
  .select('NO2_column_number_density')
  .filterDate('2019-01-01', '2019-01-30');

var tmseries = ui.Chart.image.series({
  imageCollection: data,
  region: geometry,
  reducer: ee.Reducer.mean()
});
print(tmseries);

print(tmseries.getDataTable());


texts: Area calculation: Albers Equal Area vs Pseudo Mercator

texts:
I have a question concerning area calculations in South-East Asia for West-East range that goes from India to Eastern Malaysia and a North-South range from Northern India to Indonesia (xMin,yMin 74.2203,-5.89636 : xMax,yMax 118.927,30.1183).
I found two projections that might be suitable for this geographical region:

Albers Equal Area (Conical) for S-Asia (s. http://spatialreference.org/ref/esri/102028/)

and

Pseudo/Web Mercator (EPSG:3857, s. http://spatialreference.org/ref/sr-org/7483/)

However, when I compare the results of both projections, I get large differences in total area (I used QGIS field calculator for area calculation):

620,589.41 km^2 for Albers Equal

and

689,941.23 km^2 for Pseudo Mercator

That is a difference of ard. 11%

Can anyone help me to understand the differences of these two projections and to know which one is more appropriate for my purpose?

Thanks so much in advance


texts: Using output layer names from Select By Location

texts:
This is kind of from memory as I cannot access my specific code…  I am using 10.3.1 and this will be run inside ArcMap.

Overall goal:  Find which layers in my active dataframe intersect a specific layer in the dataframe, then for each unique layer name that intersects the specific layer, output a line of text to a text file (don't need help with that part). I need to do this 5 times, each with a different intersecting layer. The dataframe will have on average 20 layers but in very rare circumstances could have up to 90 or so layers. The list of possible layers present will almost always be the same, however if the layer is empty it will not be present in the dataframe (for context in case that helps).

So, for the SelectLayerByLocation part for the first run of the tool - I think this code will accomplish that (I’ve already imported arcpy and set the workspace etc, my intersect layers are defined elsewhere in the code but I showed it again here so as not to confuse)

mxd = arcpy.mapping.MapDocument("CURRENT")
df = mxd.ActiveDataFrame
intersectlayer = “Layer 1”
layerlist = arcpy.mapping.ListLayers(df)
for layer in layerlist:
  if layer.name != intersectlayer
    arcpy.SelectLayerByLocation_management(layer, "INTERSECT", intersectlayer)

(Can I add more layer names to the “if layer.name != intersectlayer” line to weed out some returns? Or will that not matter too much in terms of efficiency?)

I am stuck on how to take the results of this query, remove the duplicate layer names (because I’m sure intersectlayer will overlap another layer multiple times) or conversely keep only unique layer names, and then pass the unique layer names with no duplicates to the part of the code where the output text is written….  Is it possible to get the results of the query, and if a result is also in the layerlist then print the message?  Not sure how to accomplish that if it's possible...

A suggestion was to getcount but that result will always be > 0 because intersectlayer would not exist if it didn't intersect anything  (this is covered in an earlier part of the script that I didn't show here).

For additional context:  if intersectlayer intersects Marsh once or multiple times, I’ll print one message “Marsh is wet” and so forth with nested If statements. I've written this code separately and it works as expected as far as I can tell - I am not requesting help with this part of the task (yet), just providing it for context.


texts: Extracting layer from MODIS data in Google Earth Engine?

texts:
I am working with MODIS imageries in GEE and I am new to its interface. So I have been able to access the land cover product of MODIS (MCD12Q1). I have added imageries from 2003 uptil 2016. How can i extract images of any particular year e.g. 2005 or 2010 in GEE for visualiztion?. There is only one image per year in this product

Here is my code for fetching and displaying the collection:

var dataset = ee.ImageCollection('MODIS/006/MCD12Q1')
              .select('LC_Type1')
              .filterDate('2003-01-01','2016-12-31')
              .map(function(img){
                var d = ee.Date(ee.Number(img.get('system:time_start')));
                var m =ee.Number(d.get('month'));
                var y = ee.Number(d.get('year'));
                return img.set({'month':m,'year':y});
              });
//var igbpLandCover = dataset.select('LC_Type1');
var palette =[
  '1f8dff',//water
  '152106',//evergreen needleaf forests
  '225129',//evergreen broadleaf forests
  '369b47',//deciduous needleleaf forests
  '30eb5b',//deciduoud broadleaf forests
  '387242',//mixed deciduous forests
  '6a2325',//closed shrubland
  'c3aa69',//openshrubland
  'b76031',//woody savanna
  'd9903d',//savanna
  '91af40',//grasslands
  '111149',//permanenet wetlands
  'cdb33b', //croplands
  'cc0013', //urban
  '33280d', //crop/natural veg. mosaic
  'd7cdcc',//permanent snow/ice
  'f7e084', //barren/desert
  ].join(',');
Map.setCenter(100,15,6);
Map.addLayer(dataset,{'min':0,'max':17,'palette':palette},'IGBP classification');


texts: Python script for opening webpage directly in QGIS

texts:
I have a shapefile with points referring to a webpage with a photograph.

Until now I define an action to open the corresponding webpage in firefox. Eg.

I have to share my project. Because I don't know the configuration of the receiver I want to create an action in Python to open the webpage into QGIS:

This works if there is just one point I have to select.  But there are several points that hide each other (more than one picture taken from the same place).

With the 'Generic-action' multiple url's opens into firefox. On the other hand: with the 'python-action' only the first url, lying on top, opens.

I tried following:

from PyQt4.QtCore import QUrl
from PyQt4.QtWebKit import QWebView
layer = qgis.utils.iface.activeLayer()
myWV = layer.selectedFeatures()
for f in myWV:
   f = QWebView(None)
   f.load(QUrl(['foto_url']))
   f.show()

This does not work.


texts: Checking whether a table is geographic or not ...?

texts: What is the best generic method to check whether an opened table is geographic or not using PyQgis ?


texts: How to get or assign vertex id to polygon vertices and use these as vertex label using pyqgis?

texts:
I want to create parcels (polygons/polyline) maps with the vertices labeled incrementally by vertex id using pyqgis/QGIS. I normally do this by converting the polygon/polyline  to point layer and adding point id attribute then overlay these points (with id labels) to the polygon layers. This is ok for one polygon, however, this becomes tedious and problematic if there are already lots of adjacent polygons.

I was wondering if there is a better way to do this using pyqgis in QGIS?

I added a sample image for reference.

UPDATE: Here's what I tried so far..

layer = iface.activeLayer()

featureId = 0

for feature in layer.getFeatures() :

    geom = feature.geometry()
    polygon = geom.asPolygon()
    outerRing = polygon[0]
    vertexId = 0

    # get the number of vertices of the outer ring
    numVertices = len(polygon[0])

    featureId = featureId + 1

    #get ring vertices
    # the -1 excludes the last vertex
    for steps in range(numVertices-1) :
        vertexId = vertexId + 1

        print featureId, ", ",vertexId, ", ", outerRing[steps] [0], ", ", outerRing[steps] [1]

This allows me to access the vertices and assign IDs and generate the list shown below:

1 ,  1 ,  -1.25423728814 ,  0.79472693032
1 ,  2 ,  -0.56120527307 ,  1.1450094162
1 ,  3 ,  0.158192090395 ,  0.606403013183
1 ,  4 ,  0.0118042738386 ,  0.30526236198
1 ,  5 ,  -0.105461393597 ,  0.0640301318267
1 ,  6 ,  -0.983050847458 ,  -0.00376647834275
1 ,  7 ,  -0.79472693032 ,  0.425612052731
2 ,  1 ,  -0.56120527307 ,  1.1450094162
2 ,  2 ,  0.020286539162 ,  1.59141143305
2 ,  3 ,  0.870711284112 ,  0.957931776096
2 ,  4 ,  0.384949214866 ,  0.718276001517
2 ,  5 ,  0.158192090395 ,  0.606403013183
3 ,  1 ,  0.384949214866 ,  0.718276001517
3 ,  2 ,  0.567873765457 ,  0.504288791391
3 ,  3 ,  0.212961505246 ,  0.230321432632
3 ,  4 ,  0.0118042738386 ,  0.30526236198
3 ,  5 ,  0.158192090395 ,  0.606403013183

I tried to combine labeling script below:

def vertexLabelling(posX, posY, vertexId) :
    palyr = QgsPalLayerSettings()
    palyr.enabled = True
    palyr.isExpression = True
    palyr.setDataDefinedProperty(QgsPalLayerSettings.Size,True,True,"8","")

    # set the placement position of the label
    palyr.setDataDefinedProperty(QgsPalLayerSettings.PositionX,True,True, posX,"")
    palyr.setDataDefinedProperty(QgsPalLayerSettings.PositionY,True,True, posY,"")

    palyr.fieldName = vertexId
    palyr.writeToLayer(layer)

    return()

vertexLabelling(str(outerRing[steps] [0]), str(outerRing[steps] [1]), str(vertexId))

However, QGIS only show one label with the last part of the loop.

If the labeling method (PAL) is applicable in this scenario, how do I modify my code?

Perhaps an alternative for me is to use the generated list as labels without having to convert it to points.


texts: arcpy.GetParameterAsText not passing arguments to script?

texts:
I have been working on an ArcGIS Python script tool which takes the co-ordinate system from an input raster and reprojects a vector to the same co-ordinate system.

The script works fine when i hard code the addresses for the workspace and input datasets. But as soon as I try to pass those as parameters in a toolbox it no longer works.

Can anyone here point out what I am doing wrong?

import arcpy

class Toolbox(object):
    def __init__(self):
        """Define the toolbox (the name of the toolbox is the name of the
        .pyt file)."""
        self.label = "Toolbox"
        self.alias = ""

        # List of tool classes associated with this toolbox
        self.tools = [Tool]

class Tool(object):
    def __init__(self):
        """Define the tool (tool name is the name of the class)."""
        self.label = "Test_tool"
        self.description = ""
        self.canRunInBackground = False

    def getParameterInfo(self):
        """Define parameter definitions"""
        param0 = arcpy.Parameter(
            displayName="Input workspace",
            name="workspace",
            datatype="DEWorkspace",
            parameterType="Required",
            direction="Input")
        param1 = arcpy.Parameter(
            displayName="Input classified raster",
            name="input_raster",
            datatype="GPRasterLayer",
            parameterType="Required",
            direction="Input")
        param2 = arcpy.Parameter(
            displayName="Input features",
            name="input_features",
            datatype="GPFeatureLayer",
            parameterType="Required",
            direction="Input")

        params = [param0, param1, param2]

        return params

    def execute(self, parameters, messages):
        """The source code of the tool."""

        # Define some paths/variables
        outWorkspace = arcpy.GetParameterAsText(0)
        arcpy.env.workspace = outWorkspace
        input_raster = arcpy.GetParameterAsText(1)
        input_features = arcpy.GetParameterAsText(2)

        output_features = outWorkspace + "\\projected.shp"

        out_coordinate_system = arcpy.Describe(input_raster).spatialReference
        proj_grid = arcpy.Project_management(input_features, output_features, out_coordinate_system)

        return

When I run the tool it provides the following error:

Traceback (most recent call last):
  File "<string>", line 75, in execute
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\__init__.py", line 1245, in Describe
    return gp.describe(value, data_type)
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 370, in describe
    self._gp.Describe(*gp_fixargs(args, True)))
OSError: "" does not exist
 Failed to execute (Tool).

it seems to me that the input variables of the parameters are not being passed correctly.


texts:
I'm pretty new to using the console and still have a lot to learn so this might come across as a really stupid question.

I imputed this code into the console:

counter = 0

# Set the path for the output file
output_file = open(r'/Users/admin/Desktop/TestFlood.txt', 'wb')
iter_layer = QgsProject.instance().mapLayersByName('test_flood')
layers = QgsProject.instance().mapLayers().values()

for point in iter_layer[0].getFeatures():
   counter += 1
   print(counter)

   for layer in layers:
       if layer.name() != iter_layer[0].name():
           for feature in layer.getFeatures():
               if (feature.geometry().contains(point.geometry())):
                   msgout = '%s, %s \n' % (point["pcds"], layer.name())
                   unicode_message = msgout.encode('utf-8')
                   output_file.write(unicode_message)

print('Finished')
output_file.close()

I am trying to run a few thousand points through each of the layers I have created and get an output on a text file that I can use to analyze further. And added a counter to see how long it takes as when I first ran this it took around a week. However, when I run it this is what I get:

Traceback (most recent call last):
File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
exec(code, self.locals)
File "input", line 1, in 'module'
File "string", line 31, in 'module'
KeyError: 'pcds'

This is probably very simple to fix but I have searched online for a while and couldn't find anything to help.

I am using QGIS 3.4.9-Madeira and my OS is MacOS 10.14.6 Mojave.


texts: Iterating in ee.FeatureCollection.aggregate_array on Google Earth Engine?

texts:
I have a shapefile I made from a list of longitude and latitudes. I extracted the latitudes and longitudes in arrays following this link. I then try to iterate over the arrays using this:

for (var i = 0; i < lats.length; i++) {
  print('lat',lats[i]);
  print('long',longs[i]);
  var point1 = ee.Geometry.Point([lats[i],longs[i]]);
  Map.addLayer(point1, {color: 'cyan'}, 'geodesic polygon');
}

but for some reason both lats.length and lats[i] (similarly for longs) are undefined. I need to iterate over longs and lats array and mark the points as (latitude, longitude) on my map.

How can I sort this issue?

EDIT: Here is the link to my script.


texts: How to find shortest route between ANY source for ALL targets?

texts:
I have a list of starting points and a list of target points for pgRouting. I want to calculate the shortest path for every pair of [source,target]. I am using a subselect to pass an array of node id's to the pg_dijkstra() function:
SELECT seq, id1 AS node, id2 AS edge, cost, linestring
  FROM pgr_dijkstra(
    'SELECT id::integer, source::int4, target::int4, st_length(linestring)::float8 as cost FROM my_ways',

    --source nodes
    array(
        select
            target as "source_ids"
        from
            my_ways
        where
            id in (
                select gid from my_sources
            )
    ),

    --target nodes
    array(
        select
            target as "target_ids"
        from
            my_ways
        where
            id in (
                select objectid from my_targets
            )
    ),
    false,
    false
  ) as di
  JOIN atene_4377.ways_atene_4377 ways
  ON di.id2 = ways.id ;

But when I do so I get the following error message:

ERROR: function pgr_dijkstra(unknown, integer[], integer[], boolean, boolean) does not exist
Hint: No function matches the given name and argument types. You might need to add explicit type casts.

However, if I replace the arrays with single node ids, the algorithm works fine. So there is something wrong about how I pass the array to the function. Any ideas how to fix this?


texts: How to add a layer using actions in PyQGIS?

texts:
I'm learning to use the python console within QGIS.

I want to add layers like iface.actionAddOgrLayer(), open projects and more by the action-function, but without any response - only the string comes up in the console like "object at 0X....."

Functions like zoom works. Did I miss something here in the API-documentation?


texts: Speed up zonal statistics with additional calculation

texts:
I have two rasters, one is a classified image and the other the corresponding probabilities for a correct classification. Now I need to calculate zonal statistics for polygon segments in the following way:

get majority of predicted class within each segment
get all probability values within this segment ONLY for those pixels that have the majority class and calculate the median probability

Here is what I do (basic code only):

rasterize polygons based on FID to create an FID-array
read classification array
read probability array
create new fields in the polygon file that shall hold the final values

# id_array: holds all FIDs as pixel values
# classi_array: holds classification (per pixel)
# prob_array: holds probabilities (per pixel)
# lyr: Layer object from ogr.DataSource().GetLayer()
# all_ids: array of unique FIDs created with np.unique(id_array[id_array != NODATA])

# this is the slow part:
for fid in all_ids:
    segment = np.where(id_array==fid)
    majority_class = np.bincount(classi_array[segment].flatten()).argmax()
    probs = prob_array[segment]
    majority_probs = probs[np.where(classi_array[segment] == majority_class)]
    med_majority_probs = np.nanmedian(majority_probs)
    feat = lyr.GetFeature(fid)
    feat.SetField(class_column, int(majority_class))
    feat.SetField(prob_column, float(med_majority_probs))
    lyr.SetFeature(feat)
    feat.Destroy()

It works, but takes forever. The question is: how can I speed up this? Is there any way to use something like multiprocessing here (with which I've never worked before) or even something completely different?
Versions (I'm bound to them due to company-wide installations, sorry):

Python 2.7
numpy 1.15.4
gdal 1.11.3.

EDIT: Now that it has run, here is the performance for the loop:
Progress: 100%|██████████| 130833/130833 [5:43:43<00:00,  6.34it/s]

EDIT 2: I guess I could somehow utilize the answer from https://stackoverflow.com/questions/20190668/multiprocessing-a-for-loop, but just can't get it to work.


texts: ESRI Javascript Web Application Hosted on Amazon Web Services

texts:
I have just deployed my first web application to Amazon S3.  This web application is a simple web mapping app using the ESRI Javascript API and consuming a web map from ArcGIS Online.

I have been having issues with this application because of the proxy page.  According to the help, in order to use ArcGIS Online web maps, you have to set up a proxy page.  I didn't do that and when I was developing the application, most components worked without issues in Chrome and Firefox.  Lately though, I have noticed that the application, which supports editing, does not apply edits from Internet Explorer.  The console gives an error message explaining that it is trying to reach a url through the proxy page and is not succeeding because the proxy page does not exist.  As far as I know, S3 does not support this proxy page because it requires an application like PHP or ASP.Net.

Has anyone had any experience with this?  I am thinking my next step may be to try and deploy this web application using AWS Elastic Beanstalk instead.  Has anyone had any luck with that?


texts:
I have images which are already separated into JPEG compressed red, green, and blue bands.  Is there a convenient way of using libgdal to pass in encoded bands without relying upon libgdal to perform the compression?  Also, is it necessary that each band of a Tiff be of the same resolution?  Due to the use of a Bayer pattern filter our red and blue bands have half the resolution of the green.

Thanks in advance.


texts: Generating 3D contours instead of 2D lines from a raster

texts:
I have difficulties generating 3D contour lines from a raster file in QGIS. It's a tif raster of an elevation model that I exported from a professor's ArcGIS project. I want to create contour lines in order to import them into SketchUP to design some engineering biological measures. Every time I use the Raster>Extract>Contours command, I only get 2D lines. I've already checked old forum entries and tried a couple of solutions, but so far nothing's worked out.

What could the problem possibly be?


texts: How to get polyline coordinates in Leaflet

texts:
Using leafletjs, I created a polyline connecting a bunch of markers:

var line = L.polyline( latlngs, {
  weight: 5,
  color: '#8EE9FF'
})

Now, I am trying to get the start and end markers when I click on a polyline:

line.on('click',(e)=>{
  // get starting marker
  // get ending marker
})

Looking at e, I see the entire list of latlgns but I just need 2 for the line segment I clicked on. Is it not possible? Should I create multiple separate lines instead?


texts: Dynamic Text in ArcMap 10.0 and arcpy - getting user name in specific format

texts:
I need to put in a specific text into a text element.  The initials of the user who built the map.

With arcpy I am trying to set the text element with dynamic text.  If I put

<dyn type="user"/>

in the user text element it gives me, correctly, the user's user name.  I want to use the user name in a python dictionary to put the users initials instead.

dictUserElmt = { "jperez" : "JP" }

Using arcpy to read the text of my text element, textElement.text, returns (again correctly)

u'<dyn type="user"/>'.

How do I get the actual output of the dynamic text "jperez" so I can use my dictionary to have the text element set to "JP".

I'm sure my wording is confusing because I am confused on how to approach.  Thanks.


texts: Manually correcting raster DEM (OpenSource)?

texts:
Is it possible, to correct a raster DEM by changeing the raster-cell-values manually, like by "painting"?

My DEM (ASTER) have one small error where the elevation is about 300m below the actual elevation. It should be a hill, and not a hole...

SRTM 4.1 is good, but has other errors (the hill is very steep and small, so SRTM "ignores" it (<90m).


texts: Overview for next/previous page in Atlas

texts:
I build a map-book for a hiking trail. On every single page is a part of the trail so you can go step by step from page to page - nothing news. But usually there should be a overview on every page for the overlaping of the next and the previous pages.

Any ideas? Until now I was working by adding the map bounding from every single Atlas page to the canvas and displaying it back in the Atlas. That is not even elegant nor is it comfortable.


texts: Changing Map Selection drawing priority in QGIS

texts:
I'm using QGIS 3.16.3
I've noticed that the drawing priority of a selection, is lower than other records in the same table.
That is - my record selection gets covered or hidden by unselected records.
I haven't noticed this before and maybe its a default setting with 3.16.3.
Does anyone know if there is a setting to change the drawing priority of map selection items?


texts:
I am trying to calculate mean water pixels within buffer regions of 10 kms using the JRC Monthly Water History data.
When I used the ee.Reducer.mean() option, I received the message: Error: Error in map(ID=00000000000000000713): Dictionary: Unrecognized arguments to Dictionary constructor.
Please let me know how to create a dictionary correctly in the function. (I am new to GEE).
var FeatureCollection = fc.map(function(feature){

  var geometry = feature.buffer(10000).geometry();

  var counts = water.iterate(function(current, previous){
    var hist = ee.Dictionary(current.reduceRegion(ee.Reducer.mean(), geometry, 30).get('water')).values();
    var count = hist.get(3);
    var date= ee.Date(current.get('system:time_start')).format('MM-YYYY');
    var values = ee.List(hist.get(0)).map(function(i){return date.cat('_').cat(ee.Number(i).format())});
    return ee.List(previous).add(ee.Dictionary.fromLists(values, count));
  }, ee.List([]));
  return feature.set("stats",counts);
});


texts:
I am using QGIS 3.6.2-Noosa. Last time I merged rasters (a couple of weeks ago) I had options on the Raster Menu that are mysteriously gone now (I would merge via Raster>Miscellaneous>Merge). Now when I click on the Raster menu I have only two options: Raster Calculator... and Align Rasters...

I am not aware of having updated QGIS and cannot find documentation of Merge on QGIS' website. Any idea what is going on?


texts: How can I simplify a 'pixelated' raster in QGIS?

texts:
I generate a Vegetation Index (VI) with a script and write the resulting 2d array into a .TIF file. When I load it into QGIS, I get the following result:

As you can see, the raster looks very 'granular'. I would like to simplify the raster by grouping large groups of pixels and removing very small 'islands' of pixels. Is there a tool in QGIS available to do this?
EDIT:
I am using the suggested Majority Filter of SAGA. In the log, I get the following error however:
QGIS version: 3.16.4-Hannover
Qt version: 5.12.8
GDAL version: 3.2.1
GEOS version: 3.9.0-CAPI-1.16.2
PROJ version: Rel. 7.2.1, January 1st, 2021
SAGA version 7.3.0 is not officially supported - algorithms may encounter issues
Processing algorithm…
Algorithm 'Majority filter' starting…
Input parameters:
{ 'INPUT' : '/home/sytze/Code/DeepLabV3Plus_implementation-single-class/my_mask.tif', 'MODE' : 0, 'RADIUS' : 1, 'RESULT' : 'TEMPORARY_OUTPUT', 'THRESHOLD' : 0 }

grid_filter "Majority Filter" -INPUT "/tmp/processing_DbxJIV/824e2db5d3434b009d4140e50b49825a/mymask.sgrd" -MODE 0 -RADIUS 1 -THRESHOLD 0 -RESULT "/tmp/processing_DbxJIV/a99c9a0147ce45679fcee5414f3d266f/RESULT.sdat"
____________________________

##### ## ##### ##
### ### ## ###
### # ## ## #### # ##
### ##### ## # #####
##### # ## ##### # ##
____________________________

SAGA Version: 7.3.0

Library: Filter
Category: Grid
File: /usr/lib/x86_64-linux-gnu/saga/libgrid_filter.so
Description:
Tools for the manipulation of gridded data.

Tools:
[0] Simple Filter
[1] Gaussian Filter
[2] Laplacian Filter
[3] Multi Direction Lee Filter
[4] User Defined Filter
[5] Filter Clumps
[6] Majority/Minority Filter
[7] DTM Filter (slope-based)
[8] Morphological Filter
[9] Rank Filter
[10] Mesh Denoise
[11] Resampling Filter
[12] Geodesic Morphological Reconstruction
[13] Binary Erosion-Reconstruction
[14] Connectivity Analysis
[15] Sieve Classes
[16] Wombling (Edge Detection)
[17] Wombling for Multiple Features (Edge Detection)
[18] Simple Filter (Restricted to Polygons)

Tool Chains: Grid Filters
Description:
Filters for grids.

Tools:
[SieveAndClump] Sieve and Clump
[grid_notch_filter] Notch Filter for Grids
[grid_list_simple_filter] Simple Filter for Multiple Grids

Error: select a tool
type -h or --help for further information

Execution completed in 0.57 seconds
Results:
{'RESULT': '/tmp/processing_DbxJIV/a99c9a0147ce45679fcee5414f3d266f/RESULT.sdat'}

Loading resulting layers
The following layers were not correctly generated.
• /tmp/processing_DbxJIV/27317bdae6e44f098ae6f732f727637c/RESULT.sdat
• /tmp/processing_DbxJIV/a99c9a0147ce45679fcee5414f3d266f/RESULT.sdat
You can check the 'Log Messages Panel' in QGIS main window to find more information about the execution of the algorithm.


texts: How to trace SQL queries sent by geo Server to Oracle database?

texts:
When I open GML view of WFS layer in geoserver, it throws this exception:

"java.io.IOExceptionORA-00918: column ambiguously defined"
I am unable to trace SQL queries sent by geo Server to Oracle database.


texts:
I'm still new into the GIS World, and here's what I'm trying to achieve. I have Yemen administrative shapefile with list of governments, and districts, I also have an Excel file with a list of migrants information. One of the fields in that Excel file is the government where that migrant was encountered.
How do I show how many migrants were encountered in each government in the map?
I could join the Excel data with the shapefile, but I didn't know how to count how many migrants with the same government and project it on the map using the symbology.


texts: Is it possible to add rings using other polygons?

texts:
I would like to "Add Rings" but using other polygons. Example:
https://www.dropbox.com/s/ue1k3b646ivqedo/holes.JPG

I would like to delete the two selected polygons and create holes under them in the big polygon where they are.

Is that possible?

Thanks!


texts: Unable to change band rendering style of rasters in QGIS 2.10.1

texts:
I am new to QGIS.

I have multiple DEM raster files (.tif) added to my map. However, I am unable to alter any of the style options under the "Band Rendering" section.

I have attempted to make changes and then save it as a new style, just clicking okay, just clicking apply, etc. Yet, every time, the original style choices are still there when I reopen the raster's properties. I am able to make changes to all other style options.

I'm particularly interesting in changing the color gradient from "Black to white" to "White to black".


texts: Using SVG icons controlled over variable order attribute field

texts:
I have about 200 different point categories. I want to order all these categories with one background SVG and one for the front. One solution is to categorize the points depending on one attribute field and after that change to rule based ordering the icons. But this caused many clicks. Is it possible to do this without so many clicks? Path and the icon name is also in a attribute field.

I don't see the function override. Now I try it, but it doesn't work. I get question marks on each position.

 Like this. What is wrong?


texts: arcpy license: import arcserver is ignored at 10.3.1

texts:
PS C:\> taskkill /F /IM python.exe
SUCCESS: The process "python.exe" with PID 14896 has been terminated.
PS C:\> python
Python 2.7.8 (default, Jun 30 2014, 16:08:48) [MSC v.1500 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> import arcserver
>>> import arcpy
>>> arcpy.ProductInfo()
u'ArcEditor'

Not sure what's up here.  I've had ArcGIS Server and ArcGIS Desktop installed simultaneously on my machine for a long time (10.1.1, 10.2.2, and now 10.3.1).  The above always worked fine to make sure arcpy scripts were using the ArcGIS Server license instead of grabbing a desktop license.  After upgrading both to 10.3.1, arcpy seems to completely ignore the server license and always look for a Desktop license.

Edited to show that no other module or process had already imported arcpy.


texts:
Is there a simple (low res is fine) polygonal data set that divides the entirety of Earth's surface into continents, with no overlap and no area excluded?

Specifically: North America, South America, Africa, Europe, Asia, Oceania and Antarctica.

I am working with a data set where most of the points are located in seas and oceans. An approximate sense of which continent all these points 'belong' to is needed.

All continent shapefiles I've seen so far trace out the coast lines, however this excludes these marine locations.

Even using country EEZ lines to map to countries and then to continents will leave much of the ocean's surface unclaimed by any continent.

I understand that such a data set might not provide a valid definition of the continents, but in my case I require all points on Earth to be associated with a continent.

Alternatively, if I must create such a data set myself, can anyone recommend a tool for doing so?

EDIT lynxlynxlynx made the interesting suggestion of using tectonic plate data:

However, for my purposes I would have to subdivide the Eurasian plate, and the ocean to the immediate west of North America would end up in the Pacific plate which is best matched with Oceania.

EDIT 2

I have a data set of lat/lng points of scuba diving locations. I wish to collate all of the dive sites at which one of my users has dived and produce a list of all the continents in which that user would reasonably think that they have dived.  As one can dive anywhere on the surface of Earth, I would ideally like to cover the earth with such a diagram.

I'm undecided about outliers such as Hawaii, which might be considered Oceania by location yet North America politically. I would tend to favour geographical location in such cases, but if a solution is easier the other way around, then that's fine too.


texts: Is there any tool or node/way tag in OpenStreetMap to obtain a POI rate as with the Google Places service? (here nodes/places have a label/tag with the rating. e.g. <rating>3.2</rating>)


texts:
I use data from NOAA for some analysis in R and I want to transform coordinates to EPSG:54004 or something useful.

There is something what they write about their coordinates.

...
POS: 29-34
 GEOPHYSICAL-POINT-OBSERVATION latitude coordinate
 The latitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where southern
 hemisphere is negative.
 MIN: -90000 MAX: +90000
 UNITS: Angular Degrees
 SCALING FACTOR: 1000
 DOM: A general domain comprised of the numeric characters (0-9), a plus
 sign (+), and a minus sign (-).
 +99999 = Missing

POS: 35-41
 GEOPHYSICAL-POINT-OBSERVATION longitude coordinate
 The longitude coordinate of a GEOPHYSICAL-POINT-OBSERVATION where values west from
 000000 to 179999 are signed negative.
 MIN: -179999 MAX: +180000 UNITS: Angular Degrees
 SCALING FACTOR: 1000
 DOM: A general domain comprised of the numeric characters (0-9), a plus
 sign (+), and a minus sign (-).
 +999999 = Missing
...

And my problem is, that I can't transform this coodrinates correctly. I use this R command:

    sc <- cbind(st$LAT, st$LON)
    ptransform(sc/180*pi,
            '+proj=latlong +ellps=sphere',
            '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

But coordinates are somehow wrong as you can see. The points should be in borders because they are meteostations from this countries. .


texts: Best way to enforce valid table names when writing to a File Geodatabase

texts:
I have an FME Workbench which reads MapInfo files and writes them to a File Geodatabase.  Generally, the workbench works pretty well, however, I've just come across an error which I though would have been handled by FME automatically.

One of the files has parentheses in the filename: AsBuilt_Crossings(Start).TAB

When FME processes this file, I get the following error:

An error occurred while attempting to create the feature class AsBuilt_Crossings(Start)_point'.

The error number from ArcObjects is: '-2147220654'.

The error message from ArcObjects is: {The table name is invalid.}

I'm guessing that the geodatabase doesn't like having the ( and ) characters.  This is something that I would have thought FME would take care of automatically, but clearly not.  Is there a way I can remove these characters (and any other invalid characters) from the name before it writes it?

I'm thinking a StringReplacer transformer with some sort of Regular Expression to replace the invalid characters with an underscore.  However, I'm not sure if I should replace the characters on fme_basename or fme_dataset.  I'm also not sure of the regular expression to use.

Am I on the right track here?


texts: Using Trace tool with BMP file in ArcGIS Desktop?

texts:
I am digitizing a map
using ArcGIS 10.

When I use trace is there any way to trace black lines from a BMP file?


texts:
I would like to have an idea about the classification of point clouds loaded in a LAScatalog with lidR. With a las file, I can do this way:
library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")
las <- readLAS(LASfile)
table(las$Classification)

How can I perform such queries with LAScatalog?


texts: Getting bathymetric slope data on ArcGIS Desktop from ASCII file?

texts:
I have a project in which I have to project XYZ ASCII file data to get a DEM. It's my first maritime map, any ideas which tools should I use?

I have already tried the Spatial Analyst Tools and haven't got any results.


texts: Calculating azimuth values along line for its segments using PostGIS

texts:
How can I calculate the azimuths along a line on its segments ?

Is there a function in PostGIS for this?

I want to calculate the min, max, sum and avg value of the azimuths for each segment.


texts: How to Stop ArcGIS API For JavaScript Drawing After Click on Map

texts:
Can you please take a look at this demo and let me know how I can stop drawing after adding the first point to the map?

as you can see I tried using the toolbar.finishDrawing();

but this is not stopping the drawing and adding new points on each click on mao

  $("#btnPoint").on("click", function () {
      toolbar = new Draw(map);
      toolbar.activate(Draw.POINT);
      map.hideZoomSlider();
      toolbar.on("draw-end", addToMap);

      function addToMap(evt) {
          var graphic = new Graphic(evt.geometry, config.symbolPointFlag);
          pGraphicsLayer.add(graphic);
      }
      toolbar.finishDrawing();
 });


texts: Select by attribute crashes ArcMap when used with search cursor, but not by itself

texts:
I have a script that I want to: run through polygons in a shapefile representing service areas > select each area > then select points from a different layer that fall within that selected polygon. (This last part isn't in my code yet.)
If I run each individual part on its own things run fine, but when I run the whole script it runs the first iteration and then I get:

Runtime Error!
Program: C:\ProgramFiles(x86)\ArcGis\Desktop10.3\bin\ArcMap.exe
R6025 - pure virtual function call"

and then it crashes.
This is my code:
import arcpy

arcpy.env.workspace = "X:\GISData\blah.gdb"

# loops through areas
with arcpy.da.SearchCursor("ServiceAreas", "Area") as cursor:
    for row in cursor:
        # select an area in the areas layer
        where_clause = "Area= '{0}'".format(row[0])
        print(where_clause)
        arcpy.SelectLayerByAttribute_management("ServiceAreas","NEW_SELECTION", where_clause)

If I just run the 'with/as' and for/in part of the code, it runs fine, and prints out a list of areas as expected.
with arcpy.da.SearchCursor("ServiceAreas", "Area") as cursor:
    for row in cursor:
        # select an area in the SCADA areas layer
        where_clause = "Area= '{0}'".format(row[0])
        print(where_clause)

If I just run the 'select layer by attribute'  line (hard-coding my where_clause), it runs fine, selecting the appropriate polygon.
arcpy.SelectLayerByAttribute_management("ServiceAreas","NEW_SELECTION", "Area= 'LR2'")

Clearing my selection between loops doesn't help. This has me stumped. Can anyone point me in the right direction? I'm open to changing the structure of my code entirely if I'm going about this in a weird way.


texts: Transit Directions: Preferred Vehicle Type?

texts:
I would like to be able to get transit directions through the JS api, that is limited to a single vehicle type (e.g. Bus). This functionality is available on maps.google.com, but how do i submit the vehicle type as part of the DirectionsRequest?.
Using the code below will provide routes via Bus, Train

google.maps.DirectionsTravelMode.TRANSIT

But if i want to restrict only to Train, how do i do it?


texts: Adding raster data to ArcGIS Online web scene?

texts:
I would like to add a number of raster datasets to an ArcGIS Online web scene, including a LiDAR DEM, a hillshade derived from it, and some aerial photography. ArcGIS Online does not support uploading any of the most common raster formats (e.g., tif), but does support rasters formatted as kml. Using the Layer to KML tool in ArcMap, I converted several of the raster datasets. They upload to ArcGIS online just fine and can be viewed or added as layers in the Map Viewer (2D). However, kml is not supported as a format in the 3D Scene Viewer as pointed out by Russ below.

I have also tried publishing the layers as Tile Layer Services without much luck. The rasters appear to upload fine, are recognized as data sources and preview at the correct extent and location in Web Scenes, but completely crash ArcGIS online and never load if added as a layer. They are in WGS 1984, which should be supported in global scenes and the files are not very large.

The source files (raster layers in file geodatabase) load and display correctly in normal offline ArcGIS Pro scenes and support use as elevation surfaces. So I'm not sure where the disconnect is when adding them to the online platform.

How do I add raster data to an ArcGIS online Web Scene?

I'm open to using any format or upload method that has worked for others, especially if they enable use as elevation surfaces in Web Scenes.


texts:
As a total Python newbie, I'm trying to install the Planet CLI  on my Windows 7 (64-bit) machine so that I can ultimately use the Planet API to clip and download imagery to my particular area of interest. When I visit the Planet Help page I can see that Planet recommends installing with the following command:

pip install planet

Into what interface do I type this command? e.g. is there a specific program I need to have downloaded first?


texts: How to use the Raster Calculator in PyQGIS (Console) in QGIS 3?

texts:
I want to use a script to do some raster calculations. I found this and this posts which seem to answer this question, however, it seems that they work for QGIS 2 and not for QGIS 3. I think the issues in QGIS 3 arise from using this line

qgis.utils.iface.legendInterface().layers()

which to my understanding gets each one of the layers in the project (each layer is defined as lyr in the loop here). In the second reference the layer is defined as bohLayer and i assume it was obtained in a similar way.

I tried using the path to a layer but obtained the following error

TypeError: str cannot be converted to qgis._core.QgsRasterLayer in this context

So my specific question is how can i call a layer in QGIS 3 to include it in a raster calculation (i.e. how to define lyr here or bohLayer here)?

Alternatively, a more general question is how to use the raster calculator from the console in QGIS 3?


texts: GRASS within QGIS - can I access grass modules through shell?

texts:
I have QGIS/GRASS installation and I would like to do my process workflow through shell in Ubuntu. It is raster vectorisation -> selection by value -> saving as shp. I have done first step by gdal command in shell. I have copied command from gui tool where it was shown. But I do not know how to access GRASS module v.extract from shell. Is there a specific command how to access these GRASS modules in QGIS processing tool through shell ?

I would like to access more than just this modul so I would like to know if it is possible or I have to use python scripts or I have to install GRASS standalone to use modules through shell.


texts: Vector Tiles: Considering a rule, change each building color of Buildings layer

texts:
I'm consuming a Mapbox hosted map with Mapbox GL (but if needed, I can use leaflet or whatever) that has the standard Mapbox street layers including OSM Buildings layer plus a custom layer uploaded by me with the shapes of the zip codes of a country.

Question:
I have 10 zip codes and I want to have each of them with a different color, so 10 colors. The buildings inside each zip code will receive their respective zip code color, so I will have 10 groups of buildings. How can I do that?

I know that I could create my own layer of buildings already classifying the buildings. But I think this won't be possible because actually, I want to use buildings of a whole country, so I would end up with a huge layer. If I manage to define the color of each building intersecting on-the-fly as they load, then I can still use OSM Buildings layer provided by Mapbox.

Any ideas?


texts: Why is QGIS reading in my CSV numeric data as string data, and how do I change it?

texts: QGIS keeps reading in my numeric data from .csv files as string data. This means I have to  figure out how to read it in with the proper data type, or convert it back to numeric inside of QGIS, neither of which I know how to do. Why is it doing this, and how do I fix it?


texts:
I have a polygon and also some point data
I want to clip the point data so that I am to see all the points that aren't in  I want to see all points that aren't in the polygon boundary

Is there a simple way to do this in QGIS ?
I don't know whether I can invert the selection .
So my aim is to get a shapefile of all points that aren't contained in the polygon itself


texts: Georeferencing Himawari-8 images with Rasterio

texts:
I was pulling some Himawari-8 images from RAMMB here and wanted to geo-reference the high res 5500 x 5500 images. I found a related question and answer using GDAL here which does work fine but I was hoping someone more familiar with GDAL and Rasterio could show me how to do the equivalent in Rasterio so I can write it fully in Python (I want to convert a lot of these images).
Here's an example of the working GDAL commands I wish to translate to Rasterio (slightly adapted from the linked question):
gdal_translate -a_srs "+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=140.7 +no_defs" -a_ullr -5500000 5500000 5500000 -5500000 full_disk_ahi_true_color_20210113014000.jpg temp.tif
gdalwarp -overwrite -t_srs "+proj=latlong +ellps=WGS84 +pm=140.7" -wo SOURCE_EXTRA=100 temp.tif Himawari8.tif

I'm hoping to get a Python snippet using Rasterio that generates the same geotif.


texts: Using Python new line (\n) with ArcGIS geoprocessing progress dialog?

texts:
I am currently using in my ArcGIS 10.3 toolbox script;

# Get the spatial reference
spatialRef = arcpy.Describe(template).spatialReference.name

#Loop through shapfiles in folder and reproject
for fc in fcList:
fcspatialRef = arcpy.Describe(fc).spatialReference.name
if fcspatialRef != spatialRef:
    arcpy.Project_management(fc, outFolder + "\\" + fc, template)

# Print shapefile Project results
arcpy.AddMessage("\n" + fc)

# Return any errors
except:
arcpy.AddMessage(arcpy.GetMessages())

Is there a way to seperate these results in the geoprocessing dialog window by a line spacing because it looks cluttered?

arcpy.AddMessage("\n" + fc)

kind of works but gets;

Running script Project...

abc.shp

xyz.shp

123.shp

How can I make the list like this please?;

![enter image description here][1]


texts:
I want to collapse a leaflet control on my map.  I have the current code:

    var legend = L.control({position: 'bottomright', collapsed:true});

    legend.onAdd = function (map) {
    var div = L.DomUtil.create('div', 'info legend');
    div.innerHTML = "
DISCLAIMER:
text

";
        return div;
    };
    legend.addTo(map);

But I still cannot get it to collapse or appear on hover like my layers selection does.  How do I need to adjust my code?  I have visited the Leaflet examples but it seems like the collapse feature is only for layer elements and not simply for a control.


texts: Modeling data in PostGIS when Several Join queries with Geometry are Needed

texts:
Suppose you have a database where you have several tables with geometric information and you have several cases in your system where you need to query the tables with a lot of joins (ie: LEFT JOIN table2 ON st_intersects(table1.the_geom,table2.the_geom)).

I believe this type of query is slow compared to a "normal" int field join, even if you use spatial indexes on the fields.

In a system where this situation arrives a lot, how do you model the data?
For example, I have many situations where you would have to query this, joining with 10 other tables.
Each table holds a different type of information, but I need to see if a row in table1 intersects with the geometry of all the other tables.

I thought about creating a many-to-many relationship, relating the ids to make the queries run faster. But then I would need a trigger to each time a new row is inserted, check which ones it intersects on each possible relationship (10, in the example) and save it. This would make the inserts a lot slower too.

Is there a common pattern to solve this type of problem on PostGIS?


texts:
I'm using eCognition developer trial, and when developing a ruleset in the process tree sometimes I make an error. However, just deleting that line of code (after extracting) doesn't undo that line in my image.

Is there a way to 'undo' only my most recent line of code?

I have only been able to delete the code then re-run my entire process tree, which takes a really long time.


texts: WMTS Trying to get list of OS tiles based on extent using OS API

texts:
I have a project and trying to implement the OS API WMTS, previously we used the OpenSpace approach but now we need to migrate.

In OpenSpace there was a method call GetMap but this does not exist in the OS API v1, however GetTile does but I need to specify the row and column for the tile along with the Zoom level.

The problem is I do not know what tile I need since I am using an extent.

Does anyone know a way of working out how I might be able to achieve this or offer an alternative?


texts: Difference between UTM Zone 15N and UTM Zone 15T?

texts:
I'm working on a project where the deliverable is a spreadsheet of point locations and needs to include Eastings and Northings in UTM Zone 15T.  I usually use EPSG:26915 NAD83 / UTM zone 15N in this study area, so I was confused about the 15T.  I did a little bit of reading about UTM Zones and see that N just means Northern Hemisphere, and that T is a division of the zone above 40 degrees latitude.

My question:  If I use EPSG:26915 when calculating the Eastings and Northings will my values be correct for the specified UTM Zone 15T coordinate system? I'm assuming this is correct but would love to have confirmation!


texts:
I have pansharpened my spot 6 image within grass 72 using i.fusion.hpf and it work well!

i.fusion.hpf pan=SPOT6_P_20151113@PERMANENT msx=SPOT6_MS_20151113.1@PERMANENT,SPOT6_MS_20151113.2@PERMANENT,SPOT6_MS_20151113.3@PERMANENT,SPOT6_MS_20151113.4@PERMANENT suffix=_hpf
i.colors.enhance red=SPOT6_MS_20151113.1._hpf@PERMANENT green=SPOT6_MS_20151113.2._hpf@PERMANENT blue=SPOT6_MS_20151113.3._hpf@PERMANENT

Now I would like to export this work in a multi-band raster. For that I have made a group.

i.group group=SPOT6_MS_20151113_hpf input=SPOT6_MS_20151113.1._hpf@PERMANENT,SPOT6_MS_20151113.2._hpf@PERMANENT,SPOT6_MS_20151113.3._hpf@PERMANENT,SPOT6_MS_20151113.4._hpf@PERMANENT

And tried to export this group... by when I open the exported file in Qgis ...

 r.out.gdal --overwrite input=SPOT6_MS_20151113_hpf output="/home/delaye/Documents/futurSahel/IMAGERIE SAT/pansharp/SPOT6_MS_20151113_hpf.tiff" format=GTiff type=Float64

the resolution is not what I'm waiting for

I missed something, but where ?


texts: Is it possible to use the symbology of an already published map service, so as to create a renderer in javascript api?

texts:
I have created a map service with fixed symbology.
I want to create a graphics layer which will use the same symbology, through a renderer. In other words, I want to use the already existing symbology which is contained in an another map layer (feature layer).


texts:
I am using Make Table View to copy data from one table to a table view object and in the process hide and/or rename some fields. When I output the data to a dbf it hides and/or renames the fields accordingly. When I output the data to a geodatabase, it hides the specified fields but does not rename the ones that should be renamed. Any ideas why this is happening?

After some further testing, it looks like it is changing the field name, and not the field alias in the geodatabase. I tried changing field.name to field.aliasName, but that did not work.

import arcpy
from arcpy import env

# Set the enviroment
env.workspace = '//GIS-SDE/ParcelData.mdb'
env.qualifiedFieldNames = False

# Get the fields from the G_AccountInfo table in ParcelData.mdb
fields = arcpy.ListFields("G_AccountInfo")

# Create a fieldinfo object
fieldinfo = arcpy.FieldInfo()

# Iterate through the fields and set them to fieldinfo
for field in fields:
    if field.name == "ACCOUNT_NUM":
        fieldinfo.addField(field.name, "ACCOUNT", "VISIBLE", "")
    elif field.name == "OWNER_NAME1":
      fieldinfo.addField(field.name, "OWNERNAME", "VISIBLE", "")
    elif field.name == "OWNER_ADDRESS_LINE1":
      fieldinfo.addField(field.name, "OWNERADDRESS1", "VISIBLE", "")
    else:
        fieldinfo.addField(field.name, field.name, "HIDDEN", "")

arcpy.MakeTableView_management ("G_AccountInfo", "AccountInfo", '//GIS-SDE/ParcelData.mdb', fieldinfo)

# Give the output table a unique name
output = '//GIS-SDE/ParcelData.mdb/cleaned_data'

# Add the table back to the gdb
arcpy.CopyRows_management("AccountInfo", output)


texts: How do I assign a variable to a SQL Server ArcSDE database's feature class?

texts: I'm using ArcGIS Desktop 10 and have no experience using Python with a SQL Server ArcSDE database.  When using a shp file or file gdb it is very simple to reference a given feature class by simply typing the file path to that file or gdb.  However, this seems to be much different when working with mdf files with SQL Server.  I can't find any help with this in the  resource center.  I did read some about making a arcsde connection but I had no luck with that.  Any help with the Python code to do this would be a great help.


texts: comparison of atmospheric correction methods

texts: I have noticed that dos methods for atmospheric corrections are in the section of radiometric correction i.*.toar on GRASS GIS. Does that means that dark object subtraction (DOS) is a method for doing radiometric and atmospheric corrections at once? Input images for that method have to be in DN format I would suspect. I want to compare the results with 6s algorithm.


texts: GeoServer WCS performance affected due to full-database query

texts:
We host a very large image set using GeoServer's imagemosaic. Our dataset uses two dimensions (time and elevation). All WMS queries work like a charm in less than 200ms (including round trip time) but WCS queries take from 5 to 11 seconds.

Looking into GeoServer's logs I've found that every time a WCS request is processed there are two calls to the database querying for the full list of values of each dimension. This is an example:

2019-03-21 20:57:02,688 DEBUG [geotools.jdbc] - CREATE CONNECTION
2019-03-21 20:57:02,689 DEBUG [geotools.filter] - exporting PropertyName
2019-03-21 20:57:02,689 DEBUG [geotools.jdbc] - SELECT distinct("date") FROM "public"."data"
2019-03-21 20:57:04,496 DEBUG [geotools.jdbc] - CLOSE CONNECTION
2019-03-21 20:57:04,505 DEBUG [geotools.jdbc] - CREATE CONNECTION
2019-03-21 20:57:04,506 DEBUG [geotools.filter] - exporting PropertyName
2019-03-21 20:57:04,506 DEBUG [geotools.jdbc] - SELECT distinct("asset") FROM "public"."data"
2019-03-21 20:57:06,604 DEBUG [geotools.jdbc] - CLOSE CONNECTION

As you might see above, that full table queries take 4 seconds in a best case scenario.

We are making the WCS request using the POST method with this XML for example:

<GetCoverage service="WCS" version="2.0.1">
    <CoverageId>test__images</CoverageId>

    <Format>image/tiff</Format>

    <DimensionTrim>
        <Dimension>Lat</Dimension>
        <TrimLow>-34.677094</TrimLow>
        <TrimHigh>-34.667741</TrimHigh>
    </DimensionTrim>

    <DimensionTrim>
        <Dimension>Long</Dimension>
        <TrimLow>-61.546615</TrimLow>
        <TrimHigh>-61.526359</TrimHigh>
    </DimensionTrim>

    <DimensionSlice>
      <Dimension>elevation</Dimension>
      <SlicePoint>203</SlicePoint>
    </DimensionSlice>

    <DimensionSlice>
      <Dimension>time</Dimension>
      <SlicePoint>2019-03-18T00:00:00Z</SlicePoint>
    </DimensionSlice>

</GetCoverage>

We tried using the GET method also (using subset parameters) but we are observing the same results.

Is there any way we can avoid WCS to make that full table query?


texts: How to enable LZW compression in a QGIS python code?

texts:
I have this simple code that batch sieves all the raster in a given directory. I found outthat the resulting output files (241MB )are pretty much larger than the input(~31MB). In ArcGIS, there is a way to enable compression like LZW, how to do it it QGIS python? Is there a way to insert it in this code?

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

input = "/data/brent/GAUP/sen2agri_two_years/"
output_dir = "/data/brent/GAUP/sen2agri_two_years/sieved//"
os.chdir(input)

for fname in glob.glob("S2*.TIF"):
    print fname
    name= list(os.path.splitext(fname))
    name[0] = name[0] + "_sieved"
    newName = "".join(name)
    print "processing " + newName
    subprocess.call(["gdal_sieve.py", "-st", "20", "-4", fname, "-of", "GTiff", output_dir+newName])


texts: Simple Fill and Line Pattern Fill in Categorized Symbols in the same layer

texts:
I have a unique layer with different polygon features. The symbology for each one is a simple fill with a random color, except if the name of the feature starts with "UR". In this case, the color is still random but the filling must be a Line Pattern.
I have to do it for QGIS using Python but I have to admit that I am a starter in programming.
I know the simple fill can be done by QgsFillSymbol.createSimple(), but how can I get the line pattern filling? I have used the QgsLinePatternFillSymbolLayer.create(), but I am sure I am missing something...
Is there a way to do it using Categorized Symbols?
Here is a little example of what I have. I know it has mistakes, but as I told you, I am just a starter and this is only the "chassis" of what I am trying to do.
if name.startswith("UR"):

    symb = QgsLinePatternFillSymbolLayer.create( { 'angle' : 30, 'distance' : 2.0, 'width' : 0.5, 'color' : color } )
    symb.setlineAngle(angle)symb.setDistance(distance)
    symb.setLineWidth(width)
    symb.setcolor(color)
else:
    symb = QgsFillSymbol.createSimple( { 'color' : color,  'outline_color' :'gray'  } )


texts:
I am trying to use the bng locator in ArcGIS Desktop. Unfortunately, I am not very successful. I follow following steps:

ArcGIS Desktop - ArcMap

Open the Catalog window
Expand the GIS Servers Folder and select Add ArcGIS Server
In the wizard, choose Use GIS Services
Click Next
Enter https://datahub.esriuk.com/arcgis/services or http://datahub.esriuk.com/arcgis/services (no authentication details are required)
Click Finish
Select the new server and choose a locator from the list.

I can add the server to ArcCatalog and I can see it in the Address Locator list. Unfortunately, after entering for example: SE338191 I get information that no address has been found.

What am I doing wrong?


texts: Is there a way to calculate a viewshed from n number of points? e.g., what is the cumulative viewshed when using each corner of a building

texts: I'm attempting to run a viewshed analysis in QGIS using r.viewshed, but the terrain is hilly and building is wide and low, so I get dramatically different results depending on which corner of the building I use (or any other part of the building).  Is there a suggested workflow for generating a viewshed which shows me the visibility across the landscape of any part of the building?  Can I generate a viewshed from multiple points?


texts: How to serve data that changes with time

texts:
This question is practically the same as an older one that concerns GeoServer; I'm interested in MapServer instead (and I don't see a satisfactory answer to the GeoServer question either).

I have some GeoTIFF files, each one being a temperature map at a specific time. Each hour I receive new data, so I create a new GeoTIFF file. I want the client to be able to ask the server for the list of times for which there is data. I can't understand how to do that.

I experiment with WCS spatio/temporal indexes, but GetCapabilities and DescribeCoverage don't seem to give me that information (I can't see the list of available data anywhere).

I'm not supposed to have a program modify the map file each time I get new data, am I?


texts: Changing area units on Line Density Feature

texts: Why can't I change the area units on Line Density Feature? Is there any solution?


texts:
I'm working with Leaflet trying to get some GeoJSON going. Everything works fine when I directly define a var as the GeoJSON I'm trying to display:

    var censustracts = {
    "type": "FeatureCollection",
    ...

That works great but fills up my script with coordinates and whatnot. The problem comes when I try to use the L.geoJson() function on the file itself. This throws an error:

    var censustracts = L.geoJson("tractsedited.geojson").addTo(map);


texts: osm2pgsql - can't find "maxspeed" tag after import

texts:
I have imported a compressed osm file ".pbf" in a PostGIS database with osm2pgsql using the following command:

osm2pgsql -c -d osm_france -U postgres -H localhost -S default.style -C 10000 france-latest.osm.pbf

Everything was fine, except that I can not find the "maxspeed" information  in any of the tables "planet_osm_roads" and "planet_osm_line." Does the PBF file contain this information? If it contains it, how can I extract it?


texts: How long should I hold onto (cache) references to IWorkspace, ITable, IFeatureClass?

texts:
After many months of programming with ArcObjects, I still feel somewhat unsure about two related issues:

For some operations, workspaces and datasets must be explicitly opened before they can be used. (For example, workspaces are opened via IWorkspaceFactory.Open(), ITables via IWorkspace.OpenTable(), etc.) On the other hand, there are operations that can be performed without having to open workspaces or datasets first (e.g. geoprocessing, or operations that accept INames as inputs).

Thus I suspect that there is a certain cost involved with opening a workspace or dataset, which can only be avoided for the latter type of operations. And because of that, I tend to keep (cache) references to the opened IWorkspace, ITable, IFeatureClass, … instances whenever I open these for the first time.

My first uncertainty is whether caching IWorkspace, ITable, etc. references could ever get me into trouble. More specifically, are there any ArcGIS operations, e.g. on tables, that require that there be no live references to "open" tables?

(I'm aware of only one potential problem source at the moment: If a workspace or dataset gets deleted, all cached references to it should be thrown away.)

My second, related uncertainty concerns how long I should cache workspace and dataset references (at most). Let's say I'm writing an ArcGIS extension; would it be OK to open required workspaces when the extension is loaded, and only reset the references when it is unloaded? Would it even be OK if I did the same for datasets?

Links to ESRI reference documentation (or insider knowledge) is most welcome.


texts: Interpolating multiple years of daily precipitation data to generate daily time series?

texts:
I am trying to use ArcMap ver. 10.2 to create several models using ModelBuilder. The model will be created to:

interpolate 7 years of point measured precipitation data (long integer format at a daily time step) using various surface-fitting methods (e.g. inverse distance weighting, polynomial, and kriging)
clip the output raster files to subbasins
do zonal stats on the clipped raster files
append the zonal stats into a table.

I have built the models, and they all run well (Fig. 1).

The models only work for one day at a time. However, I need to run each model for over 2,000 iterations (one iteration for each day).
I've tried using the iterator tools, but there is not a tool that will work with the way I have the attribute table set up (Fig. 1).
I need to generate daily time series' of interpolated precipitation data by looping my models for multiple iterations.
How can I make my models loop for multiple iterations, one iteration for each set of field values (i.e. column in attribute table)?

Additional background information follows:

The precipitation data are stored in a point file with a large attribute table. The attribute table has a field for each day precipitation was observed. The fields are labeled in sequence "PPT". Each field has multiple point measured daily precipitation values (long integer format), one point measurement for each location (n=10) where precipitation was measured.

Fig. 1. Model schematic and attribute table used to interpolate precipitation data.


texts:
I am trying to get labels to follow the linestring. I have inserted the following sld text into a geoserver style, but the label is displayed at the horizontally at linestring centroid

          <sld:TextSymbolizer>

       <sld:Geometry>
          <ogc:Function name="centroid">
           <ogc:PropertyName>shape</ogc:PropertyName>
          </ogc:Function>
       </sld:Geometry>

        <sld:Label>
          <ogc:PropertyName>road_name</ogc:PropertyName>
        </sld:Label>

        <sld:Font>
          <sld:CssParameter name="font-family">MS Shell Dlg 2</sld:CssParameter>
          <sld:CssParameter name="font-size">6.0</sld:CssParameter>
          <sld:CssParameter name="font-style">normal</sld:CssParameter>
          <sld:CssParameter name="font-weight">normal</sld:CssParameter>
        </sld:Font>

        <sld:LabelPlacement>
          <sld:LinePlacement>
            <sld:PerpendicularOffset>10.0</sld:PerpendicularOffset>
          </sld:LinePlacement>
        </sld:LabelPlacement>

        <Halo>
            <Radius>
                <ogc:Literal>2</ogc:Literal>
            </Radius>
            <Fill>
                <CssParameter name="fill">#FFFFFF</CssParameter>
                <CssParameter name="fill-opacity">0.85</CssParameter>
            </Fill>
        </Halo>

        <sld:Fill>
          <sld:CssParameter name="fill">#000000</sld:CssParameter>
        </sld:Fill>

       <sld:VendorOption name="followLine">true</sld:VendorOption>

      </sld:TextSymbolizer>

I thought that by adding the VendorOption -> "followLine" the label would be displayed inline with the line. Is there something that i have missed in the documentation to make this work?