texts: ModelBuilder Collect Values output values to a table?

texts:
I am working with Build Balanced Zones tool in ArcGIS Pro, which uses a genetic algorithm to find an optimal solution to a spatial problem (such as creating zones with equal population). It uses a parameter called generations (the more generations you use, the better the solution will be) to generate those zones. There isn't much information about what is the optimal amount of generations because it largely depends on the problem. By using the default settings, it shows that it found the optimal solution, but by running those same parameters again, the solution is different. I believe it gets stuck in local maxima. I found that by increasing the generations to a fairly high number (which makes the tools running time fairly time-consuming) it finds the same solution every time. I want to find out how many generations is it reasonable to use in order to find the solution which is the same every time I run it.

I believe there are many ways to find the solution, but right now I created a model with ArcGIS ModelBuilder, which is missing a couple of steps in the end.

It starts with a For loop, which increases the number of generations used in the iteration. It outputs the generated solution.

From there, the solution is dissolved based on the generated Zone_ID in order to compare it against the perfect solution which was previously created using a very high number of generations. Ideally, it should compare it against the next generation, but it is difficult to use two iterators together. The Spatial Join used is Are Identical To, which outputs the join count of 1 if the polygons are geometrically identical and 0 if they are not.

Then I use the Summary Statistics tool which finds the sum of the join count field. Generating a table with one row in the process. I use a Get Field Value tool to find that value.

I am calculating that value based on an expression. Right now, I am using smaller sample data and I know that if that value equals 4 then it is the right generation. I could use a Stop tool right there, which ends the iteration, but I found that the result could also be false positive. Genetic Algorithms use randomness and it could have been just luck that the two solutions match perfectly.

Hence, I need to collect the values and apply some sort of pattern finder script to it (if the expression is true ten times in a row for example). Currently, I have collected those values with Collect Values, but it does not let me output it to some sort of a table, for which I could apply further expressions. The values are just there, and I can't seem to do anything with them. I can clearly see that the expression became true (value of 1) in a visible pattern as the generation number started to increase. Is there any way I can create a view or table with those values in the same iteration?


texts: Creating multiple points, using Python AddIn?

texts:
import arcpy
import pythonaddins
class point(object):
    """Implementation for point_addin.tool (Tool)"""
def __init__(self):
    self.enabled = True
    self.shape = "NONE"
def onMouseDownMap(self, x, y, button, shift):
    pt = arcpy.Point()
    pt.X = float(x)
    pt.Y = float(y)
    ptGeoms=[arcpy.PointGeometry(pt)]
    arcpy.CopyFeatures_management(ptGeoms, r"C:\test.shp")

There is a problem in my code. While i use this add-in in arcgis, it creates a new point on where i click, but the problem is, it deletes last point i created. I can not save my last point.


texts: Stripes/Lines appear on QGIS rasters in different zoom levels, why?

texts:
The stripes rotate as zoom level changes

Is there an option to stop this effect?


texts: Unable to save to in GeoJSON format with the correct projection

texts:
I have a shapefile (Polygon) in QGIS that I wish to save in GeoJSON in order to display it with the mapbox JavaScript API. Original projection is in POSGAR 97.

My issue is that it does not reproject when I set projection to WGS84. Even if in the T.O.C I set the projection to WGS84 and Save As. Even if I uncheck the on the fly re-projection.

This is how coordinate looks like:

"coordinates": [ [ [ 324345.0, -3579435.0 ]

It should looks like:

"coordinates": [ [ [ -32.43450, -35.794350 ]

Curious think is if I Save As in KML, it works and the point of the polygon looks like:

<coordinates>-64.866399861691903,-32.338025763844477...


texts: Google Earth Engine exported image with missing records in some bands

texts:
I've been exporting an ImageCollection (Sentinel-2 data) on Google Earth Engine by stacking the images with the toBands() function. Before, during image collection I've used a function to mask clouds using the Sentinel-2 QA band (maskS2clouds).
Download of the stack (consisting of a single band) worked without problems but, when processing data outside GEE (some simple Grid Statistics for Polygons on SAGA) I noticed that some bands had no values at all.
Has anyone ever experienced something similar?
If so, how did he resolve the problem?


texts: Network analysis tool (service area from layer) doesn't seem to work for path type fastest

texts:
I am trying to map the areas which are a 15 minute drive from recycling centers.
I am having issues when using the Network analysis service area (from layer) tool when using the fastest path type. No matter what I input into the cost field my entire road network is mapped out, whether I say 60 seconds or 900 seconds (15 mins). I have tried a number of different times and they all give this same result.
I tried using the shortest path type and this seemed to work fine.
I have used this tool before in another project without issue so I am not really sure what the problem is.


texts: How to integrate QGIS 3.4 with the DE theme

texts:
I just got the upgrade to QGIS 3.4 from the ubuntugis-unstable PPA. One thing that immediately stands in this version is that it ignores the DE theme, instead yielding an unhealthy white-based theme. In the figure below you can compare it with GEdit, an application that automatically integrates with the DE.

I have been using QGIS for more than a decade and this is the first version that does not integrate with the DE. I have gone through the Settings menus and could not find any option to activate integration with the DE. How can it be done?

Update: as requested by Erik Friesen, where is the output of qt5ct:

$ qt5ct
Configuration path: "/home/user/.config/qt5ct"
Shared QSS paths: ("/home/user/.local/share/qt5ct/qss", "/usr/share/i3/qt5ct/qss", "/usr/local/share/qt5ct/qss", "/usr/share/qt5ct/qss", "/var/lib/snapd/desktop/qt5ct/qss")
Shared color scheme paths: ("/home/user/.local/share/qt5ct/colors", "/usr/share/i3/qt5ct/colors", "/usr/local/share/qt5ct/colors", "/usr/share/qt5ct/colors", "/var/lib/snapd/desktop/qt5ct/colors")

And the main settings in qt-config:


texts:
I am giving this for intersection of a polyline using Leaflet draw plugin with polygons feature class. Here is my code
function intersectPolyByPolyFC(poly,fcPoly){
        var fgp=[];//variable for feature group
        var bbPoly=turf.bboxPolygon(turf.bbox(poly));//polygon representation of the bounding box of a reference polygon
        for(var i=0;i<fcPoly.features.length;i++){//loop through the feature collection of polygons
            var bb=turf.bboxPolygon(turf.bbox(fcPoly.features[i]));//variable bb represents each polygon for feature collection
            if(turf.intersect(bbPoly,bb)){//condition to check if the two bounding boxes intersect
                var int= turf.intersect(poly,fcPoly.features[i]);
                if(int){
                    int.properties=fcPoly.features[i].properties;
                    fgp.push(int);
                }
            }
        }
        return turf.featureCollection(fgp);
    }

**Function Definition**
var bufLine= turf.buffer(line,5,{units:"kilometers"});
                        var intBUOWL = intersectPolyByPolyFC(bufLine,
                        lyrBUOWL.toGeoJSON());
                        L.geoJSON(intBUOWL,{style:{color:'black', width:5}}).addTo(mymap);

I am not getting any error on console and moreover my intersection is not performed on the map. What is could the problem be with my code?


texts: Storing shapefiles in recently created folder using ArcPy?

texts:
I'm trying to create a folder and export the following shape files to the new folder. I got it to create the folder but I can't seem to get the shapefiles to the new folder.

import arcpy
import datetime
from arcpy import env

env.workspace = "C:/UPDM.gdb"
out_folder_path = "C:/output"
nowstart = datetime.datetime.now()
YearMonthDay = nowstart.strftime("%Y_%m_%d")
out_name = "WebTMS_BU_" + YearMonthDay + ""

arcpy.CreateFolder_management(out_folder_path, out_name)

arcpy.FeatureClassToShapefile_conversion(["P_Leaks"],
                                         out_name)


texts: Why don't my featureLayer labels display when using timeSlider widget (Javascript API for ArcGIS 4.14)

texts:
You only see all the labels when the slider has gone through the full time extent; at initialization there is nothing, they then start rendering as the slider plays.

Anyone know what the issue is? See my codepen here: https://codepen.io/cvricella2/pen/eYNxdQg


texts: QgsMapCanvas layers disappear on canvas panning/zooming

texts:
I am trying to write a plugin that will show an image based on its footprint in the vector layer.
Currently, I have a problem with QgsMapCanvas. Layer in QgsMapCanvas disappear when I pan/zoom canvas. Here's my code:
myPlugin.py
def startOblique(self):
        self.dlgSec.canvas_dlg.show()

        self.dlgSec.canvas_dlg.enableAntiAliasing(True)

        newRaster = QgsRasterLayer(r'D:\basic\FRONT\B002_0891.tif','Raster1')
        print(newRaster)
        print(self.dlgSec.canvas_dlg)

        self.dlgSec.canvas_dlg.setLayers([newRaster])
        self.dlgSec.canvas_dlg.setExtent(newRaster.extent())

        self.dlgSec.canvas_dlg.refresh()
        self.dlgSec.canvas_dlg.waitWhileRendering()

def run(self):
        if self.first_start == True:
            self.first_start = False
            self.dlg = test1Dialog()
            self.dlgSec = test2Dialog()

        project = QgsProject.instance()
        # show the dialogs
        self.dlg.show()
        self.dlgSec.show()

        self.dlg.oblique_starter.clicked.connect(self.startOblique)

When I click oblique_starter button I get this:

So, everything is fine by now. Raster is on the canvas.
When I try to pan/zoom on the canvas raster disappears so I get this:

Extra info:
When I execute simple code to show canvas and raster layer inside python console on QGIS I can pan and zoom layer and it doesn't disappear:
canvas = QgsMapCanvas()
canvas.enableAntiAliasing(True)

newRaster = QgsRasterLayer(r'D:\basic\FRONT\B002_0891.tif','Raster1')

canvas.setExtent(newRaster.extent())
canvas.setLayers([newRaster])

canvas.resize(QSize(1366, 1064))
canvas.show()

What can I do to prevent the layer from disappearing from the canvas in my plugin?


texts:
I'm trying to retrieve the metadata of an items on my portal.

I've been trying that : /sharing/content/items//info/metadata/metadata.xml but it's not working.

I have Portal for Arcgis 10.3 and I've installed geoportal server.


texts: QgsNetworkReplyContent.content() empty after HTTP error

texts:
I have a custom Flask backend and a QGIS plugin as client. I make a request with the QgsNetworkAccessmanager with blockingPost() to an endpoint which returns a JSON response on either HTTP status.
That works fine for HTTP 200, but when an API error occurs, the content() of the QgsNetworkReplyContent is just empty (b''). I'm fairly sure I'm missing smth on the Flask side, because in other plugins requesting other (non-Flask) APIs I can correctly see the content error message.. I tried diving into the source code to see what the Flask response is missing, but it's quite hard since the whole network module mixes Qt and Qgs components quite heavily.
EDIT: This is the (quite standard) error handler in Flask:
from werkzeug.exceptions import BadRequest

@api.errorhandler(BadRequest)
def handle_user_bad_post_error(e: BadRequest):
    return {'error': e.description}, 400

EDIT2: forgot to say that all other requesting methods I tried (curl, Postman) did response with the correct error message


texts: Python scaling issues and possible subprocess conflicts?

texts:
I have always found such great analytical and programming expertise on StackExchange in the past that I hope you can help me.

Background: I have a project that involves spatially joining and interpolating traffic volume data across a road network on a continental scale. In order to handle the data processing, I have broken the data down by state and then into 10km processing tiles. States will run on parallel processes, which call a subprocess module to actually do the data processing in a 10km area at a time until the state level script is finished. The interpolation process starts with a road segment that has joined traffic data and searches connecting street segments with matching road name until it finds another segment with data. Then the code interpolates for segments in the middle.

The problem: When I run a single 10km tile or even ~20-30 tiles (for the same state), I get interpolation values that are expected given the algorithm I defined. When I start processing the entire state (thousands of tiles) with the very same code, I will get interpolation values for the same road that exceed the endpoint counts by orders of magnitude in either the positive or negative direction (and how could my traffic volume go negative, right?). Considering that the calculations will change for the same road depending on how much data is being processed, I feel that this is somehow more of a scaling or instability issue with Python, and not a bug in the code logic.

I was initially doing all the looping by referencing nested ArcGIS searchcursors (agh, experience backs up all the forum recommendations to NOT do that). So, I used an ArcGIS 10.0 searchcursor just once to read the street segment geometry and needed attributes into a Python dictionary for faster "lookups" during looping. The performance over ArcGIS searchcursors (even the new arcpy.da module) is amazing. I can now process the entire continent's worth of street data in two days. But now I'm wondering if Python dictionary tables can become instable if they get too large, or if the subprocess is not properly releasing memory objects before returning to the master script. There isn't any data passed from the subprocess back to the master; it's written to a permanent geodatabase table before closing. I wondered initially if the problem was running parallel processes, but I get the same unstable behavior for a single process on the state-level. I'm really not an expert Python programmer, though I've been effective in adapting other people's solutions to large data processing. This feels like one of those questions that a real developer could answer, as opposed to a code "pirate" like me.

I don't know that posting all my scripts will be helpful, as I think I'm looking for someone to point me in the right direction to research possible issues/solutions with python/subprocess instability, dynamic changing of values in dictionaries, etc. I think if I had a good starting place on issues that cause that kind of behavior, I could problem solve from there (hopefully).


texts: Reading an ESRI XML Workspace Document with python or other FOSS

texts:
I've recently run into data in "ESRI XML Workspace Format," I have never run into this format before, and can't find any reference to using it with FOSS tools, or any documentation on the format the geometry is stored in.

My goal is to read this data with python so i can processes it with Shapely, but i'm open to any FOSS solution for reading it.

Possibilities:

GDAL/ORG
QGIS
A python library

I am not interested in any solution using ArcPY.

Here is a sample of one of the records from the file.

<Record xsi:type='esri:Record'>
    <Values xsi:type='esri:ArrayOfValue'>
        <Value xsi:type='xs:int'>830</Value>
        <Value xsi:type='xs:string'>Bankhead Jones</Value>
        <Value xsi:type='xs:string'>Bureau of Land Management</Value>
        <Value xsi:type='xs:int'>2</Value><Value xsi:type='xs:short'>2</Value>
        <Value xsi:type='esri:PolygonB'>
            <Bytes>MwAAIPjj2b8URlrAYG9mt8AORkAMJpONhkNawCCTb/2OEUZAAQAAAB0AAAAAAAAADCaTjYZDWsBI
rXx4ERFGQARv5XiHQ1rAoBiidZoQRkBsfDzg2UNawJhVg8OaEEZAdH/3udpDWsBQwafHIxBGQExt
zpPbQ1rAWMCLy6wPRkB0XYht3ENawAgssM81D0ZAfFQUxi5EWsDwTfbjNw9GQLwrdowvRFrAYG9m
t8AORkAoYYfggURawGhch8PCDkZAzJQIGtREWsCAvSJixA5GQDQkkkEeRVrAcEtS1sUORkBgPvVP
HUVawBBlbkw8D0ZAxLtjfhxFWsC4m6sWsw9GQEyZnXBvRVrAcBxRpLIPRkDYdtdiwkVawICY2y+y
D0ZAQA/We8FFWsDoPX9bKRBGQPjj2b8URlrAAEunFSkQRkB0XYD8E0ZawGCUfmegEEZAfLMHORNG
WsCgpGm4FxFGQMQuldoSRlrAIJNv/Y4RRkAcx2WVv0VawDCgmNuOEUZAfGzcrGtFWsCokEu5jhFG
QNgRU8QXRVrA2HcdlY4RRkCobz1jz0RawDBQxWyNEUZAhJ/aDn1EWsDIuTIajBFGQMDNg519RFrA
AAfqihURRkAELQNSK0RawADsDaUTEUZAkI5lBtlDWsCoD2G/ERFGQAwmk42GQ1rASK18eBERRkAA
AAAA</Bytes>
        </Value>
        <Value xsi:type='xs:string'>{9E4D50CD-773F-426A-9D85-ACEB53A6651E}</Value>
        <Value xsi:type='xs:string'>BFAHRER</Value>
        <Value xsi:type='xs:dateTime'>2015-06-23T20:07:05</Value>
        <Value xsi:type='xs:string'>BFAHRER</Value>
        <Value xsi:type='xs:dateTime'>2015-06-23T20:07:05</Value>
        <Value xsi:type='xs:double'>0.00059271317454369463</Value>
        <Value xsi:type='xs:double'>0.12326579683440851</Value>
    </Values>
</Record>

The data I am trying to read is from http://www.blm.gov/wy/st/en/resources/public_room/gis/datagis/state/state-own.html


texts: How much precipitation must fall within a catchment area to fill up a sink’s entire volume assuming that the runoff conditions are perfect

texts:
I am performing an exercise and I need to answer these question, any hints would be extremely useful.

How can it be calculated how much precipitation that must fall within a catchment area to fill up a sink’s entire volume assuming that the runoff conditions are perfect (i.e. the entire runoff from the precipitation reaches a local sink)?


texts:
I am trying this example with different input files. My aim is to merge the country and point tables into one table using sjoin from geopandas and write to a CSV file. I took a look at 1.
When I run the second line of code below, It takes really long time.

My Input point file is about 62kb and my country file is about 8kb is that the reason for the time consumption ?
Is there any other efficient way to write to a CSV? Or is there any mistake in my code below ?

merged_gdf = gpd.sjoin(point_gdf, country_gdf, how="right", op="within")
merged_gdf.to_csv("plswork.csv",index=False,mode='w')


texts: Obtain new polygon outline from 2 intersected polygons?

texts: There are 2 polygons intersected. I can get a intersected small polygon by using "ST_INTERSECTION(geom, geom)". But I want to outer join the polygons like the uploaded image.


texts: Identify polygons with multiple attributes

texts:
I'm trying to identify polygons that contain both of two different land cover types from the same landcover layer. I cannot figure it out. Ive tried combinations of select by attribute and select by location, to no avail. I've also tried to intersect and query based on a single layer, but still no luck. At this point i'm sure it's something fairly simple and i'm going to kick myself once i know the answer.

I'm using ArcGis 10.4.1


texts: accessing main window in QGIS when processing modeller is active

texts:
While working on a model in the processing modeler I can't enter the main GUI.

Every time I want to check some features or layer properties I have to leave the modeller, do the check and reopen the model I'm working on.

Is there a way to swith between the modeller and the main window?


texts: Export to CAD - preparing the DWG seed file

texts:
I'm preparing a seed file in order to export to CAD and maintain all the attributes at all fileds. after creating the block definition
I've export to CAD (Right click) using the seed file I made (after adding RefName field {string 255} and insert the block name) but the attributes didn't inserted to their Autocad fields.

Here are my 3 Questions -

1.How to connect the GIS fields to Autocad blocks attributes ? is the RefName field is enough ?
2. Is there an automatically way to create seed file with blocks from a feature class ?
Is it possible to import CAD block from exiting feature class instead of manually type all the fields  ?
3. What the difference between Autocad Civil to Autocad Map in terms of working with GIS layers ?


texts:
I made a viewpoint symbol in QGIS by styling a point to be be a circle and adding some lines that are rotated at different angles, and it ended up like so

In a perfect world I would rather this symbol would be a polygon-style (rather than a point-style), so it would be the correct size at any zoom (without creating zoom-dependant symbols).

My problem is that I would like a half-circle, like this:

and a more-than-half-less-than-full circle like this:

Any suggestions on how I can create a style like this in QGIS?


texts: Error for point.in.poly function in the spatialEco package in R stating arguments imply differing number of rows


texts:
Possible Duplicate:
  How to add GeoServer wfs layer on OpenLayers?

I'm really new to this and am making a map for a university project, I can't get my WFS layer to show up on my map, the WMS layer shows up fine, just no WFS layer, here is my code:

Ext.onReady(function() {

var map = new OpenLayers.Map('map');

map.addControl(new OpenLayers.Control.OverviewMap());
map.addControl(new OpenLayers.Control.LayerSwitcher());
map.addControl(new OpenLayers.Control.PanZoomBar());
map.addControl(new OpenLayers.Control.MouseToolbar());
map.addControl(new OpenLayers.Control.KeyboardDefaults());

var base = new OpenLayers.Layer.WMS(
"UK base map",
"http://localhost:8080/geoserver/wms",
{layers: "topp:GBR_adm2"}
);
map.addLayer(base);

var renderer = OpenLayers.Util.getParameters(window.location.href).renderer;
renderer = (renderer) ? [renderer] : OpenLayers.Layer.Vector.prototype.renderers;

var wfs = new OpenLayers.Layer.Vector("WFS", {
strategies: [new OpenLayers.Strategy.Fixed()],
protocol: new OpenLayers.Protocol.WFS({
version: "1.1.0",
url: "http://localhost:8080/geoserver/wfs",
featurePrefix: "topp",
featureType: "reptiles_amphibians_locations",
featureNS: "http://www.openplans.org/topp",
srsName: "EPSG:4326",
geometryName: "the_geom"
})
});
map.addLayer(wfs);

var mappanel = new GeoExt.MapPanel({
region:"center",
height: 400,
width: 600,
map: map,
title: 'UK Distribution'
});

var layerList = new GeoExt.tree.LayerContainer({
text: 'Species',
layerStore: mappanel.layers,
leaf: false,
expanded: true
});

var layerTree = new Ext.tree.TreePanel({
title: 'Map Layers',
maxWidth: 500,
region: "east",
collapsible: true,
collapsed: true,
root: layerList
});

mainPanel = new Ext.Panel({
height: 600,
renderTo: "mainpanel",
layout: "border",
items: [mappanel,layerTree]
});

});

I keep reading that it could be due to the fact i'm not using a proxy, but its on localhost:8080, so does it still need one? It could be a geoserver problem. Its been driving me mad for weeks, any help would be much appreciated.

Thanks

Will


texts: Multi-temporal stack does not include all bands

texts:
I made a stack of images for a multitemporal series for 7 months (one per month with 5 bands) with gdal, but the stack that comes out to me consists of a single band for each month, probably the first.
How can I include all 5 bands for the 7 months?

Here's my code:

outtif = 'C:\Users\tack.tif'

tifs = [
     r'C:\Users\19_4_2018.tif',
     r'C:\Users\4_5_2018.tif',
     r'C:\Users\23_6_2018.tif',
     r'C:\Users\18_7_2018.tif',
     r'C:\Users\22_8_2018.tif',
     r'C:\Users\21_9_2018.tif',
     r'C:\Users\31_10_2018.tif'
   ]

outds = gdal.BuildVRT(tifs, separate=True)
outds = gdal.Translate(outtif, outds)
outds=None


texts: Add buffer only on one specific side of polygons

texts:
I have a shapefile with multipart geometry (multipolygons), namely communities (blue on the picture below). There are also lakes (white), but for these I don't have a geometry, they're just "empty".

I want to conduct a point-in-polygon analysis with another point dataset. Each point should be assigned to a community or be deleted when it falls outside a community (e.g. in a lake).
The problem is that there are many points for lake-bordering communities that are very close to the shore. Because of the strong generalization, these points get accidentally deleted even though they would actually belong to a community.

The only thing I can think about is somehow to extend the shore line of lake-bordering communities by about 100m, but I'm wondering whether there is some way of doing this automatically in QGIS. I have no clue though I would approach this problem.


texts: Input feature mask name with variable - error 865

texts:
I have 73 raster files, with 73 corresponding shapefiles. I have to extract each of these raster by their corresponding shapefiles. I tried to call the shapefile using %Name%, since it has the same name, but it gives me an error.
My shapefiles are from Thiessen Polygon, and in the end I have to build a mosaic with all these pieces of raster.
How do I solve this problem?


texts: NAD83 to NAD83(CSRS) transformation for Northwest Territories, Canada

texts: I am looking to determine the correct Grid Shift Binary (gsb) file for transforming data from the NAD83 original to NAD83(CSRS) geodetic datum for the NWT in Canada. I was able to find reference for a NAD83-98.gsb file for this purpose for Quebec, and I believe I saw another for Ontario, and other provinces, but nothing for NWT.


texts:
I'm using GeoServer for serving map in the GeoJSON format. I din't find any working example.  Here is my code snippet.

<!doctype html>

  <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.0.1/css/ol.css"
    type="text/css">

  <title>OpenLayers example</title>

URL is working perfectly in the postman but not in the OpenLayers. Did I miss anything?


texts: Activating new Snapping feature in ArcGIS Engine while editing features and vertices?

texts:
I'm working on new Snapping Capability of ArcGIS Engine.
I want to use snapping capability using ISnappingEnvironment class while users are editing a feature or moving vertices.
I tried IPointSnapper.CacheShapes() and set the FeatureSelection of IMap as argument of type GeometryBag  but I dont know when I should cache features and clear it or update cache.
We are using ControlsEditingEditToolClas control for editing features and I want to use new snapping with this tool or other available options.


texts: Geocode free APi for internal application

texts:
I need opinions / experienced use on this issue:

Which API can be used to geocode an internal application of a company that meets these requirements:

Can be used for company internal application unlimitedly
Returns the latitude and longitude of the address more precisely and not atrealada the vertices (as with OSM)

Right now I'm researching on Yahoo Maps on license terms. I want opinions from someone who has passed this need. My application is to design routes using OSRM and Leaflet. I'm currently using Nominatim Geocode which returns me the lat / long in accordance with the vertices when the address is in the middle of the street - this causes problems.


texts: How to modify features in a WFS transaction with the Geoserver scripting extension

texts:
I'm trying to modify a feature in a WFS insert transaction using the Geoserver Scripting Extension with python. Basically I'm able to change the feature attributes, but I'm unsure of how to pass this feature along the pipeline so that the modified attributes will end up in the database.

def preInsert(inserted, req, context):
    iterator = inserted.features()
    while iterator.hasNext():
        feature = iterator.next()
        feature.setAttribute('userid', 'placeholder')

The above code doesn't have the userid field populated once it gets to the database.

Really I'm just confused as to whether or not I need to return something from these hook functions or if I'm allowed to mutate the objects themselves because the documentation is sparse and examples are few and far between.


texts: PyQgis question: How to change or set active layer?

texts:
Developing a QGIS plugin, I found that if no layer is selected in the Layers panel I can iterate over all layers and do anything to them, but "activate" them.

Could someone give me a hint please?


texts:
I am using arcpy.CopyFeatures_management(pointGeometry, "output") to create singles points.

Is it possible to create those single points on an existing output without deleting the output before the overwrite?

The method description says:

  The feature class which will be created and to which the features will
  be copied. If the output feature class already exists and the
  geoprocessing overwrite option is set to true, the output will be
  deleted first.


texts: Standalone error using Python 3: layer.isValid is False

texts:
I am trying to create a standalone application using Python 3.7 and QGIS3.6. on Windows 10. I have installed QGIS using the OSGEO4W method, in the 64-bit flavour.

My issue is that when I try to load a layer, it always says that the layer is not valid, as determined by layer.isValid() = False.

current paths are:

  pythonpath=C:\OSGeo4W64\apps\qgis\python
  path=C:\OSGeo4W64\apps\qgis\bin;C:\OSGEO4~1\apps\Python37;C:\OSGEO4~1\apps\Python37\Scripts;C:\OSGEO4~1\apps\Python27\Scripts;C:\OSGEO4~1\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\system32\WBem

I have read several answers on this, but none of them have worked for my so far; mostly the answers have been to do with changing setPrefixPath, which I have done in several different ways to no avail. I have also changed the pythonpath, but that did not work either.

What else can I try?

My script:

import os, sys, time, csv
os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'c:\osgeo4w64\apps\Qt5\plugins'
os.environ['PATH'] += r';c:\osgeo4w64\apps\qgis\bin;c:\osgeo4w64\apps\Qt5\plugins'

sys.path.extend([r'c:\osgeo4w64\apps\qgis\python',r'c:\osgeo4w64\apps\Python37\Lib\site-packages'])

from qgis.core import *
from qgis.gui import *
from qgis.utils import *
from qgis.core import QgsProject
from PyQt5.QtCore import *
from PyQt5.QtGui import *

#   Instantiate the QGIS Application
app = QgsApplication([], False)   # second argument to False disables the GUI

#   Update prefix path
app.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
app.initQgis()
#print (QgsApplication.showSettings())

#   File to study area (vector)
vector_source = "path/to/inputs StudyArea.shp"

#   Read layer from file
layer = QgsVectorLayer(vector_source, "Vector_Layer", "ogr")

# Test if the file path points to a valid file
print ("File Path points to a file: ", os.path.isfile(vector_source))
#  Test if the layer is valid
print ("Output of layer.isValid(): ", layer.isValid())

I finally found the answer by changing by setting the os.environ['QGIS_PREFIX_PATH']. See @Matt Brauer’s answer to PyQGIS QgsVectorLayer() Loads Invalid Layer in Standalone Python 3 Script


texts: Using Leaflet to consume google tiles without labels

texts:
Using Leaflet, it is possible to consume google maps tiles with the code below:

    var googleStreets = L.tileLayer('http://{s}.google.com/vt/lyrs=m&x={x}&y={y}&z={z}&',{
        maxZoom: 20,
        subdomains:['mt0','mt1','mt2','mt3']
    });

Is it possible to consume the same service but without labels? Does google provide a tiling service for this purpose?

I can easily do that with cartodb basemaps because they provide different services for this reason:

    var basemapCartoNoLabels = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_nolabels/{z}/{x}/{y}.png', {
        attribution: '&copy;OpenStreetMap, &copy;CartoDB'
    }).addTo(map);

    var basemapCartoOnlyLabels = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_only_labels/{z}/{x}/{y}.png', {
        attribution: '&copy;OpenStreetMap, &copy;CartoDB',
        pane: 'labels'
    }).addTo(map);

My goal is to show the following layers:
- google streets
- my geometry json layer
- labels on top of everything


texts:
I'm trying to convert polygons to points with ArcGIS 10. My polygon layer has 2 attributes that I want to give to my point layer. I'm using the Feature to Point function from the Arctoolbox, but it always results in certain polygons having two or more points. I tried without and with the inside option. Anyone have a solution?

Thanks


texts: Create a layer selector when making a map at runtime for CartoDB

texts:
Using the same theory to create a search box at runtime, (such as here), I wrote this code:

var l = cdb.vis.Overlay.create('layer_selector', map.viz, {})
l.show();
$('#map').append(l.render().el);

But this doesn't work as anticipated. The layer selector doesn't fill with any layers. It tells me I have 0 layers.

Can someone tell me how to get a layer selector at runtime, using the cartodb layer selector?

EDIT:
I've taken the bits of my code related to this:

function main() {
var map = new L.Map('map', {
    zoomControl: false,
    center: [40.7142, -74.0064],
    zoom: 11
});
L.tileLayer('http://tile.stamen.com/toner/{z}/{x}/{y}.png', {
    attribution: 'Stamen'
}).addTo(map);

cartodb.createLayer(map, {
    type: 'cartodb',
    options: {
        user_name: '',
        table_name: '',
        extra_params: {
            map_key: ''
        }
    },
    sublayers: [{
        sql: "SELECT * FROM roads",
        cartocss: '#roads{line-color: #FF6600;line-width: 2;line-opacity: 0.7;}'
    }]
})
    .addTo(map)
    .on('done', function (vis,layer) {

    var l = cdb.vis.Overlay.create('layer_selector', map.vis, {
        layer_selector: true
    })
    l.show();
    $('#map').append(l.render().el);

}).on('error', function () {
    console.log("some error occurred");
});
}


texts: How to transfer attribute values in a One to Many relationship in QGIS?

texts:
How could one go about recreating this Arcmap workflow and result in QGIS?

https://esriaustraliatechblog.wordpress.com/2015/06/22/spatial-joins-hidden-trick-or-how-to-transfer-attribute-values-in-a-one-to-many-relationship/


texts: postgis ST_Contains with ST_Buffer performance issues

texts:
I have a neighborhoods layer and a parcels layer. I want to find how many parcels are in each layer.

Initially, I used this query:

SELECT neighborhoods.name, COUNT(parcels_kcgis.geom) AS total
FROM neighborhoods LEFT JOIN parcels_kcgis
ON ST_Contains(neighborhoods.geom, parcels_kcgis.geom)
GROUP BY neighborhoods.name

This query returns counts from all 200 neighborhoods in less than a second.

However, I noticed that it misses some parcels on the edges of the borders. I thought a good way to do this would be to buffer the neighborhood geometry slightly:

SELECT neighborhoods.name, COUNT(parcels_kcgis.geom) AS total
FROM neighborhoods LEFT JOIN parcels_kcgis
ON ST_Contains(ST_Buffer(neighborhoods.geom, 10), parcels_kcgis.geom)
GROUP BY neighborhoods.name

Unfortunately this takes approximately 30 minutes to run. Using a where clause to only get one neighborhood takes 4-5 seconds.

I am trying to use this information in a web application and need a much quicker response time. Is there a better way to do this?


texts: Where can I find a map with the location of towns and cities in Southern Sudan?

texts:
For a project I have some name of towns and cities in the Southern part of Sudan, the region South Kordofan. Many of these towns I can't find on Google Maps, but I need to know their location.

Example of towns I can't find are: Kuma, Mendi, Kalkutta, Um Serbeba and Al Nugra.

Hope someone knows a good source?


texts: How to update/refresh fields/attributes after modifying SpatiaLite/SQLite database?

texts:
I have modified the columns in a SpatiaLite database table (by creating a new table with a renamed column, populating it with the old table, deleting old and renaming the new table to the old).

I have run VACUUM on the database, re-added the layer anew, refreshed and reconnected to the database, and also tried restarting QGIS, but I am unable to get the renamed (and reordered) columns to show correctly. In both the Attributes Form (and table) and Source Fields the old columns are used. The database manager does show the current structure.

Is there some cache to clear? How can I get QGIS to pick up the changes?

Update: oh, I see there is a table geometry_columns_field_infos which isn't shown in the Database Manager, which contains a list of fields. Is there a way to rebuild it, or do I have to do this manually?

Update 2: well, I found a workaround to update the statistics:

UPDATE geometry_columns_statistics set last_verified = 0;
SELECT UpdateLayerStatistics('geometry_table_name');

This indeed updates the geometry_columns_field_infos table, but it doesn't really solve the issue. When I try to use the SpatiaLite layer in a virtual layer, I get Query preparation error on PRAGMA table_info(_tview): no such column: geometry_table_name.columnname.


texts: Inserting tags (multi-values) in field of attribute table with QGIS?

texts:
I am setting up a QGIS project using a Geopackage (or shapefile) concerning the areas named with some historical toponyms and writing a brief description of them.

In the attribute table I would like to include a field in the form of "tags" (such as Youtube videos or this forum), where I can insert more keywords referenced in a single field to assign a "kind" to the object and facilitate the management of data.

The question is that an element in my conceptual model can have more genres.

For example, an element can have keywords at the same time: "restaurant", "hotel", "garden".

In a field in the attribute table I am allowed to insert only one of these values.

How could I enter multiple values in a field as an alternative to a static text system "restaurant, hotel, garden"?


texts: Random sampling of classification in Google Earth Engine

texts:
I have been working on a supervised classification and am now trying to place random points to do an accuracy assessment. To do this I am needing to convert each class in my classification to its own geometry as I want to place a certain amount of training points in each class. I have not been able to successfully convert the classification regions into geometries yet, but I assume it is a simple fix. I thought it could be done with image.geometry(), but it made a geometry over the entire world instead of just a single class.
var improved = classified_svm.updateMask(mask_improved);
Map.addLayer(improved, {}, "improved");
var improved_clip = improved.clip(bir);
var improved_geo = improved_clip.geometry(); <--- This is where it goes wrong

Link to code: https://code.earthengine.google.com/5005a440fa695d5aa286e226f8f26ff5
Link to bir asset: https://code.earthengine.google.com/?asset=users/anna_odell/bir_Boundary2005
Link to points asset: https://code.earthengine.google.com/?asset=users/anna_odell/trainingpoints3


texts: Does ArcGIS tracking server 10.0 supports rest url?

texts:
I'm using ArcGIS Tracking server 10.0 .

I know  Tracking server 10.1 has web connector but 10.0 version does not.

Does ArcGIS tracking server 10.0 supports rest url ?

Is it possible to use tracking services 10.0 version in the web or mobile device using rest or json ?

How can i do it ?


texts: Why gdalinfo don't show subdatasets in sentinel-2 imageries?

texts:
In this article Sentinel-2 products have 'Subdatasets' where I can open the 10 meters resolution bands. But in some imageries there isn't this sector.

For example in imagery 'de4388a9-13e1-4c5c-affe-515c9a2277c5' isn't subdatsets.

gdalinfo S2A_OPER_MTD_SAFL1C_PDMC_20170602T042051_R135_V20151209T075322_20151209T075322.xml

Driver: SENTINEL2/Sentinel 2
Files: S2A_OPER_MTD_SAFL1C_PDMC_20170602T042051_R135_V20151209T075322_20151209T0
75322.xml
Size is 512, 512
Coordinate System is `'
Metadata:
  CLOUD_COVERAGE_ASSESSMENT=7.5755689585225765
  DATATAKE_1_DATATAKE_SENSING_START=2015-12-09T07:53:22.031Z
  DATATAKE_1_DATATAKE_TYPE=INS-NOBS
  DATATAKE_1_ID=GS2A_20151209T075322_002420_N02.04
  DATATAKE_1_SENSING_ORBIT_DIRECTION=DESCENDING
  DATATAKE_1_SENSING_ORBIT_NUMBER=135
  DATATAKE_1_SPACECRAFT_NAME=Sentinel-2A
  DEGRADED_ANC_DATA_PERCENTAGE=0
  DEGRADED_MSI_DATA_PERCENTAGE=0
  FOOTPRINT=POLYGON((...))
  FORMAT_CORRECTNESS_FLAG=PASSED
  GENERAL_QUALITY_FLAG=FAILED
  GENERATION_TIME=2017-06-02T04:20:51.000005Z
  GEOMETRIC_QUALITY_FLAG=PASSED
  PREVIEW_GEO_INFO=Not applicable
  PREVIEW_IMAGE_URL=Not applicable
  PROCESSING_BASELINE=02.04
  PROCESSING_LEVEL=Level-1C
  PRODUCT_START_TIME=2015-12-09T07:53:22.031Z
  PRODUCT_STOP_TIME=2015-12-09T07:53:22.031Z
  PRODUCT_TYPE=S2MSI1C
  PRODUCT_URI=822_2017-06-02T03_59
  QUANTIFICATION_VALUE=10000
  RADIOMETRIC_QUALITY_FLAG=PASSED
  REFERENCE_BAND=B1
  REFLECTANCE_CONVERSION_U=1.03002388757155
  SENSOR_QUALITY_FLAG=PASSED
  SPECIAL_VALUE_NODATA=0
  SPECIAL_VALUE_SATURATED=65535
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

In this imegery( scihub.copernicus.eu/dhus/odata/v1/Products('ab904126-1e9f-44d4-a4ca-a86deb3e86a6') ) gdalinfo works fine and there is subdatasets.
Sorry about my english.

GDAL 2.2.1, released 2017/06/23


texts: .e00 files - can they be converted to CSV files and used in QGIS?

texts:
I am having trouble opening an .e00 file in QGIS. (QGIS keeps locking up, and I have to "force-quit".) The file is political administrative boundaries in China.

Would it be feasible to convert the .e00 file to a CSV and import that CSV? 'Cause the .e00 file, when I open it with a text editor, looks like a spreadsheet document.

I am using QGIS 2.12 Lyon on a Mac running OS X 10.10.5.

Edit: I do not see how my question is a duplicate of "How to read Arc/Info e00 vector files (no binaries) with Python?" [I'll call that "Question A".]
 First, Question A says "I can easily read e00 files with QGIS..." but in my case, I cannot easily read e00 files. I have only partial success displaying their data.
Second, I am asking whether I can turn my e00 file into a CSV file and then use the CSV file in QGIS. Question A makes no mention of CSV files.
Third, I do not mention Python. I have no idea how to use Python.
I have read Question A and the answers given, and I do not understand what they are talking about. For example, when I run "ogrinfo --formats" I get "-bash: ogrinfo: command not found".


texts: Remote sensing product's raster extent not agreeing with specified projection in R

texts:
I want to plot remote sensed sea ice concentration data for Antarctica, in Program R. The product website (http://cersat.ifremer.fr/data/products/catalogue) says the data are in Antarctic polar stereographic projection. However, the extent of the rasters (thus the x/y (lat/lon) coordinates) does not agree with this projection.
Preamble. This works, and shows what I expect:
This is a bathymetry map of the Antarctic Peninsula, in WGS84 projection.

I can reproject this map to Antarctic polar stereographic projection. In R:
Polar.prj = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "
Site.polar = projectRaster(from = site.WGS84, crs = Polar.prj)
I end up with a raster with extent values in the “negative millions” (longitude) and “positive millions” (latitude) (see map image, and raster details below)

> Site.polar
class      : RasterLayer
dimensions : 1321, 3217, 4249657  (nrow, ncol, ncell)
resolution : 460, 1090  (x, y)
extent     : -3026365, -1546545, 922441.7, 2362332  (xmin, xmax, ymin, ymax)
crs        : +proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

Sea ice concentration problem:
I wish to use sea ice concentration data from http://cersat.ifremer.fr/data/products/catalogue.
Data are available here: (ftp server) or here (3 files I downloaded; < 1 MB each)
This product contains Antarctic polar sterographic 12.5 km resolution grids of sea ice concentration.
Every day of the year has its own .nc file. I load 3 files (3 days of sea ice concentration data) in R as follows:
library(raster)
library(ncdf4)

setwd("Computer/covariates/seaice_IFREMER")

#first import all files in a single folder as a list
rastlist <- list.files(path = getwd(), pattern='.nc$', all.files=TRUE, full.names=FALSE)

# stack the raster data
icerasters <- stack(rastlist)
Here I get warnings that
[1] "vobjtovarid4: **** WARNING **** I was asked to get a varid for dimension named nj BUT this dimension HAS NO DIMVAR! Code will probably fail at this point"

#Plot sea ice data from a single day
plot(icerasters[[1]])

This looks strange! It is a mirror image of Antarctica!
# Use raster::flip to change it:
icerasters <- flip(icerasters, direction='y')
plot(icerasters[[1]])

Better. This looks like the Antarctic Polar Sterographic Projection.

But I still have a problem.
The 'values' of the x/y (lat/lon) axis (in other words the extent of the raster) do not correspond to values I associate with the polar stereographic projection.
The units of the polar stereographic projection are meters (just as the units for WGS84 are in degrees).
This means the projected bounds (extent) of the polar stereographic projection should be in 'millions', like those of bathymetry in the preamble. This is not what appears on the ice map.
Let us look at the raster information:
print(icerasters)
class      : RasterStack
dimensions : 664, 632, 419648, 106  (nrow, ncol, ncell, nlayers)
resolution : 1, 1  (x, y)
extent     : 0.5, 632.5, 0.5, 664.5  (xmin, xmax, ymin, ymax)
crs        : NA

I can set crs manually, since it is not recorded in the imported raster:
Polar.prj  = "+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs "
crs(icerasters) = Polar.prj

But this does not change the values of the extent / axes. I need these to correspond to those mapped in the bathymetry layer.
How do I do that?
The IFREMER sea ice data are derived from NSIDC 12.5 km data. I have exactly the same problem downloading data from NSIDC (https://nsidc.org/data/ae_si12). These data are stored as .he5 files. Whether .nc or .he5, I always get the warning:

**** WARNING **** I was asked to get a varid for dimension named nj BUT this dimension HAS NO DIMVAR! Code will probably fail at this
point"

when importing the data from this sea ice product. Perhaps this is the root of the problem.


texts: Drawing line between beginning and end points with lat/long data in ArcGIS Desktop?

texts:
I have a table 500 entries that each have a beginning latitude, beginning longitude, ending latitude and ending longitude.

I would like to draw a line between the beginning lat/long and the ending lat/long of each entry (row).

What is the simplest way to accomplish this?


texts: Python in ArcGis, cannot find what is wrong with my syntax

texts:
I am using Feature class to Feature class to batch convert some .dwg files. However python does not like the last statement (name).

It's probably something very simple, any idea where I'm going wrong?


texts: Is there an established method to declare "Forward" on a linear feature

texts:
Is there an accepted way to establish a Forward direction for a linear feature regardless of the ordering of the points used to store the feature?

An example would come from Roadway Design in the United States and other countries in which "stationing" is determined to be positive if the end point is east of the begin point.  (For ties in the easting value, I suppose we compare the northing value.)

I am looking for a generalized solution, though, which would be able to handle linear features which happen to have identical begin and end points.


texts: gdal_calc - How to calculate where raster A is not equal to raster B?

texts:
I have two binary rasters with values NaN and 254. I want to generate a new raster that consists of the pixels that have value NaN in raster A and 254 in raster B and vice versa - those should be 1, the rest should be NaN again.

I want to use gdal_calc.py in a bash script for that. My code snippet looks like this:

gdal_calc.py -A S2_2015_AOI_NDBI_bin.tif -B S2_2015_AOI_NDSI_bin.tif --outfile=S2_2015_AOI_diff_auto.tif --calc="1*(A!=B)" --NoDataValue=0

But all it returns is a new raster with NaN everywhere.

I also tried --calc="1*(not_equal(A,B), --calc="1*(logical_or(logical_and(A==254, not_equal(B, 254)), logical_and(B==254, not_equal(A, 254))))", and several other combinations I could come up with. Nothing worked yet.

Strange thing is, it works flawlessly in the QGIS Raster Calculator with the formula ("S2_2015_AOI_NDBI_bin@1" != "S2_2015_AOI_NDSI_bin@1")  = 1, but I want this to be automated. Also, this calculator throws 0 instead of NaN, which is not suitable for further processing as I plan it.

The two input rasters were previously generated by:

gdal_calc.py -A S2_2015_AOI_NDBI.tif --outfile=S2_2015_AOI_NDBI_bin.tif --calc="254*(A>-0.15)" --NoDataValue=0

and

gdal_calc.py -A S2_2015_AOI_NDSI.tif --outfile=S2_2015_AOI_NDSI_bin.tif --calc="254*(A>0.4)" --NoDataValue=0`

The rasters all have the same CRS.
I'm running GDAL v. 1.11.3, it's a bit outdated, might this be the problem?

UPDATE

I solved the problem by skipping the step of creating the "binary" input rasters with 254 and NaN values (see question's second and third code blocks). The code looks like this now:

gdal_calc.py -A S2_2015_AOI_NDBI.tif -B S2_2015_AOI_NDSI.tif --outfile=out.tif --calc="1*(logical_and(A>-0.15, B<0.4))" --NoDataValue=0

This delivers the result I wanted.

While this is a suitable workaround for me personally, it still does not answer the original question.


texts:
Possible Duplicate:
  Problem with Orfeo Toolbox integration in QGIS

I am trying to do is how is shown in this link: http://www.youtube.com/watch?v=XcAQDFtkbdk

I did that steps after steps and the system returned me this message:

Then I checked the Sextante log and I read this information:

1- OTB execution console output
'"C:\Users\Owner\sextante\tempdata\otbcli_MeanShiftSegmentation"' is not recognized as an internal or external command,
operable program or batch file.

2- SAGA execution commands
io_gdal 0 -GRIDS "C:\Users\Owner\sextante\tempdata\1352576101.2845.sgrd" -FILES "C:\Users\Owner\sextante\tempdata\modelertest139.tif"
shapes_grid "Vectorising Grid Classes" -GRID "C:\Users\Owner\sextante\tempdata\1352576101.2845.sgrd" -CLASS_ALL 1 -CLASS_ID 0 -SPLIT 1 -POLYGONS "C:\Users\Owner\sextante\tempdata\modelertest142.shp"

3- SAGA execution console output
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>set SAGA=C:\Users\Owner\sextante\tempdata\
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>set SAGA_MLB=C:\Users\Owner\sextante\tempdata\modules
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>PATH=PATH;C:\Users\Owner\sextante\tempdata\;C:\Users\Owner\sextante\tempdata\modules
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>saga_cmd io_gdal 0 -GRIDS "C:\Users\Owner\sextante\tempdata\1352576101.2845.sgrd" -FILES "C:\Users\Owner\sextante\tempdata\modelertest139.tif"
'saga_cmd' is not recognized as an internal or external command,
operable program or batch file.
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>saga_cmd shapes_grid "Vectorising Grid Classes" -GRID "C:\Users\Owner\sextante\tempdata\1352576101.2845.sgrd" -CLASS_ALL 1 -CLASS_ID 0 -SPLIT 1 -POLYGONS "C:\Users\Owner\sextante\tempdata\modelertest142.shp"
'saga_cmd' is not recognized as an internal or external command,
operable program or batch file.
C:\ProgramData\Microsoft\Windows\Start Menu\Programs\OSGeo4W\OSGeo4W>exit

In my opinion I should have problem with the configuration of the OTB in QGIS, but I do not know how to fix it.

I will show you my  Sextant options and configuration to be sure I have it properly.

OTB applications folder                  C:\Users\Owner\sextante\tempdata\
OTB command line tools folder       C:\Users\Owner\sextante\tempdata\

SAGA folder                                 C:\Users\Owner\sextante\tempdata\

Any help to resolve this issue?

Thank you in advance

Hilario


texts: Plotting points using list of MGRS coordinates on spreadsheet into QGIS

texts:
I'm very much a beginner at QGIS, but need to learn something fast.
I am trying to project a table of "MGRS" coordinates from an excel spreadsheet into QGIS. I need a quick and easy way to do this, either through multi point layers or however else it is possible.
My question is similar to this one...
Import MGRS vector (old projection) to QGIS

I know that QGIS seems to be very much in favor of X/Y or lat/Long coordinates, however I'm working strictly with MGRS, and I need QGIS to recognize the MGRS coordinates in that column of my spreadsheet.  Essentially if I add more MGRS points to that spreadsheet, I want to be able to plot them quickly in case of updates.

I don't need an Excel spreadsheet that converts MGRS coordinates in their separate sections of MGRS into lat/ long or x/y coordinates.  I've already found and downloaded protected spreadsheet for conversion, but I need to be able to plot points on the map from MGRS coordinates in the location column of my spreadsheet.

I don't just want a conversion website.  I'm not trying to find out what the point I'm looking for is.

I've already downloaded lat long tools, but that too has a coordinate converter, but I don't think that's the tool I need.
UTM coordinate import using csv table

I've tried using a delimited text layer, but it seems to only use lat long coordinates in separate boxes.  every instructional video on importing a spreadsheet talks about making a "delimited text layer". but delimited text layers take x/y coordinates as well, and I'm not sure how to implement MGRS coordinates.

I'm also not trying to create and project MGRS grid squares on my map.
http://gis.stackexchange.com/a/261144/49538

Can anyone give me a simple step by step walkthrough of how I'd do this?


texts: How do I make a user input array or list in python?

texts:
I don't know which one I need actually because I am relatively new to python. What my goal is to have a stand alone script to run multiple buffers based on user input. For example, I would like to ask the user how many buffers and then ask them what those buffers are.

How many buffers do you have? 3
What distances would you like to buffer? 50, 150, 250

I want to get all the input information before running the buffer tool. I'm guessing I need some type of loop.


texts: How do you create a buffer for a specific feature in QGIS?

texts:
Fairly competent ArcGIS user here, very new to QGIS.

What I am trying to do is create a one-mile buffer around the highlighted parcel in QGIS 2.14 (from a .shp).  Seemed straightforward enough: select polygon, click 'use only selected feature' (said polygon), input buffer distance in meters (I will reproject if needed, but I haven't even got that far yet), browse to output folder, and..."No feature selected, please uncheck 'use selected' or make a selection"

I did select a feature. It's right there, highlighted.  Or is there some other way to 'select a feature' that I don't know about?  Attached is a screen capture of what happens when I thought I had it figured out.


texts: Color a layer of the map based on the location of dots on another layer

texts:
I created a map as below:

The map is created using the following script:

ct <- readShapeSpatial("housect_37800_0000_2010_s100_census_1_shp/wgs84/housect_37800_0000_2010_s100_census_1_shp_wgs84.shp")

ct_mod <- fortify(ct,region="SLDLST10")
chart <- ggplot(data=ct_mod,aes(long,lat))
chart <- chart + scale_x_continuous(limits=c(-73.8,-71.7),breaks=seq(-74,-71,0.1))
chart <- chart + scale_y_continuous(limits=c(40.9,42.1),breaks=seq(40,43,0.1))
chart <- chart + geom_polygon(fill="grey80",aes(group=group))
chart <- chart + geom_path(color="white",aes(group=group))
# coord_equal make x-y same scale
chart <- chart + coord_equal()

chart1 <- chart + geom_point(data=hosp.list,aes(x=coord_x,y=coord_y))

This is a map of the Connecticut Census State Legislative District Map with several hospitals located, you can see a few regions with a dot and others without. How can I fill the regions with dot using one colour and regions without dot using another colour?


texts:
I have this code and it works fine:

<script>
  var raster = new ol.layer.Tile({source: new ol.source.OSM()});

  var wkt = 'LINESTRING (30 10, 10 30, 40 40)';
  var wkt2 = 'LINESTRING (0 0, 5 5, 10 10, 8.646 50.4654)';

  var format = new ol.format.WKT();

  var feature = format.readFeature(wkt, {dataProjection: 'EPSG:4326', featureProjection: 'EPSG:3857'});
  var feature2 = format.readFeature(wkt2, {dataProjection: 'EPSG:4326', featureProjection: 'EPSG:3857'});

  var vector = new ol.layer.Vector({source: new ol.source.Vector({features: [feature]})});
  var vector2 = new ol.layer.Vector({source: new ol.source.Vector({features: [feature2]})});

  var map = new ol.Map({
    layers: [raster, vector, vector2],
    target: 'map',
    view: new ol.View({
      center: [2952104.0199, -3277504.823],
      zoom: 4
    })
  });
</script>

Is there a simple way to define an own style (e.g. color, thickness) for each wtk (in code wkt and wkt2)?


texts: How to set QGIS options in standalone pyqgis application

texts:
I have a standalone application which uses the QGS api to view and query map layers. One layer is a WMS-C layer which automatically caches tiles which is great. I want to also activate a setting which is present in QGIS but I can't work out how to implement it in my app, namely:

Use render caching... (Settings -> Options -> Rendering)

I found the QSettings key for this in QGIS by browsing the output of:

from PyQt4.QtCore import QSettings
s = QSettings().allKeys()
for k in s:
    print k

The one I need to access/emulate/implement is:

Qgis/enable_render_caching

Where can I activate/implement these settings in my standalone app? Is this the setting that is read by QgsMapCanvas.setCachingEnabled(bool)?


texts:
I need to find temporal rolling means (or other stats - max, min, etc.) in stars spatial obects. For stars raster cubes, this would be directly analogous to the calc(rasterBrick, function(x) movingFun(x, 3, mean)) available in the raster package. I could just use raster instead of stars if all I had were raster cubes. However, I'm really liking the ability of stars to handle vector data as well, and I have some multipolygon data with time slices in stars objects that I would also like to get rolling statistics on. For example, I might try to calculate the 3-day rolling mean for the values in each polygon.
I can use something like aggregate(polygonStars, by = 'weeks', FUN = mean, na.rm = TRUE) to get the mean over discrete time units, but aggregate aggregates rather than rolls. I've tried using st_apply, but haven't had success with feeding it a rolling function like RcppRoll::roll_mean (I'm unclear what exactly it's doing when I feed it a rolling function since it wants a single output, but it ends up replacing a dimension instead of the attribute values).
I've tried just shoving the objects into the rolling functions from RcppRoll and zoo, but they expect vectors, and I haven't found a way to tell them to operate along the 3rd dimension. Perhaps that's possible? I suppose I could take the stars objects apart into vectors, do the rolling calculations, and then put back together, but that would be cumbersome, error prone, and slow.
Ideally, I'd remain in stars, but if there's not the possibility of getting rolling temporal means from stars objects, I'll switch to raster for the rasters, and would appreciate a suggestion for the polygons with temporal layers.
To illustrate using the example from https://r-spatial.github.io/stars/, I'd like to be able to get the rolling means for each pixel in the raster cube prec and each polygon in the vector cube a, as indicated by the last lines of this code (which don't work because I'm using functions from raster to show what I'm looking for).
# Modified from the data raster and vector time series example at
# https://r-spatial.github.io/stars/ to illustrate

library(stars)
library(tidyverse)

# Get the rasters: these look similar to what my data has
prec_file = system.file("nc/test_stageiv_xyt.nc", package = "stars")
(prec = read_ncdf(prec_file, curvilinear = c("lon", "lat"), ignore_bounds = TRUE))

# Get north carolina
sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf"), "nc.gpkg") %>%
  st_transform(st_crs(prec)) -> nc # transform from NAD27 to WGS84
nc_outline = st_union(st_geometry(nc))

# Do the aggregate to make polygons with temporal layers
a = aggregate(prec, by = nc, FUN = max)

# Rolling means over time
# These don't work because the functions are from `raster`, but this is the sort
# of thing I want to do.
rastRoll <- calc(prec, function(x) movingFun(x, 3, mean))
vectRoll <- calc(a, function(x) movingFun(x, 3, mean))


texts: Making clip feature optional in ModelBuilder?

texts: I want to make the clip feature an option, and also a model parameter so that the one who will use it will have the option to add it or skip it.


texts: QGIS 3 Python Load all layer from geopackage

texts:
I am loading all layers from a geopackage in a QGIS 3 canvas in Python using a list with all layernames and the following code:

databasepath = "X:/test.gpkg"
layer = ["1","2","3","4"]
for i in layer:
   iface.addVectorLayer(databasepath+"|layername="+i, i, 'ogr')

Is there a better way to load all layers from a geopackage beside the list as the layer names or count may change?


texts:
I recently downloaded QGIS 2.8 and opened an existing project in it to check how it handled it and make any changes required. This project includes several raster heatmap layers.

My problem is that many of these raster layers seem to have shifted by a significant degree (but not all of them so additionally confusing). At first I thought it was a styling issue but then realised they were being displayed outside of the visible area.

They all have the same CRS set (EPSG: 27700, British National Grid), the project is set to this also.

I can see no differences between the settings for the affected raster layers and the single unaffected raster layer. Some of the affected layers are .ntf and some are .tif. The one unaffected layer is a tif.

Re-importing the layers doesn't seem to help. Creating a heatmap from scratch leads to the same problem, very odd!


texts: Editing the colors of a 'Point Symbol' in ArcGIS Pro 2.3

texts:
I've been transitioning from ArcGIS Desktop 10.6 to ArcGIS Pro 2.3. I believe that I've a handle on many of the changes that I need. At least I've been able to find documentation that helps. However, I have a fiddly bit that I've been unable to solve.

I am trying to change the colors of the point symbol, 'Shield 1'. It comes as a dark Red and Blue, whereas my client would greatly prefer a pastel color scheme. When one selects the 'color' option, the main color (the bottom part of the shield) can be changed. I have not been able to find a way to change the top part of the shield, (the pointy bit). I feel that I must be overlooking the possibility, but I've been playing around with it for quite awhile.

Is this possible?

If you've any advice or a resource I've overlooked I'd love to know about it.

(I apologize for the double post. I asked two questions in my first post.)


texts: Displayng attribute text in layout Text element gives error: Error in setting text

texts:
I am new in Python scripting. I have a geodatabase with layer Test and this layer has field Name and Text element on layout named TEname. I am trying to show selected feature attribute Name value in Text element, but I could not do it. I can print the value for selected feature, but can't show it in Text element.

Here is the code I use:

import arcpy

mxd = arcpy.mapping.MapDocument("current")
TEname = arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT","Name")[0]

for row in arcpy.da.SearchCursor("Test","name")
     print row

...
(u'Land plot 1',)

for row in arcpy.da.SearchCursor("Test","name"):
    TEname.text = row

...

        TEname.text = row

  Runtime error Traceback (most recent call last): File "", line 1, in
  File "c:\program files
  (x86)\arcgis\desktop10.5\arcpy\arcpy\arcobjects_base.py", line 89, in
  _set return setattr(self._arc_object, attr_name, cval(val)) RuntimeError: TextElementObject: Error in setting text


texts: Filter parent layer with child table in 1-N relation in QGIS

texts:
I have a layer with a related (1-N) table and I need to filter the parent layer with values of the child table.
I have tried with the Query builder of the Provider feature filter in the layer properties, but I can't see any way to do that.
Is it possible to do this?


texts: Adding x,y labels in QGraphicsview using pyqt4?

texts:
Here is my sample program. When I need to draw a line for x axis and y label for y axis .so can any one help me how to add a line and labels to scene.i tried different ways but i didn't get the proper output.

How can I add labels to this scene?

Given below is my sample code:

import sys
from pyface.qt import QtGui, QtCore
# class ScanView(QtGui.QGraphicsView):
#     def __init__(self,X=5, Y=5,parent=None):
#         super(ScanView, self).__init__(parent)
class DemoApp(QtGui.QMainWindow):
    def __init__(self, parent=None):
        super(DemoApp, self).__init__()
        self.w= QtGui.QGridLayout()
        self.v= QtGui.QGraphicsView()
        self.w.addWidget(self.v)
        self.widget = QtGui.QWidget()
        self.widget.setLayout(self.w)
        self.setCentralWidget(self.widget)
        self._squares = []

        n_rows, n_cols = 3, 2
        squareLB = 50
        label = QtGui.QLabel("xaxis")
        label1 = QtGui.QLabel("yaxis")
        self._scene = QtGui.QGraphicsScene()
        mytext1 = QtGui.QGraphicsSimpleTextItem('label')
        self._scene.addItem(mytext1)
        mytext2 = QtGui.QGraphicsSimpleTextItem('label1')
        self._scene.addItem(mytext2)
        self.v.setScene(self._scene)

        self.pen = QtGui.QPen(QtCore.Qt.DotLine)
        self.pen.setColor(QtCore.Qt.red)

        width, height = (2 + 2)*squareLB, (3 + 2)*squareLB
        self._scene = QtGui.QGraphicsScene(0, 0, max(708, width), height)

        p = squareLB if width > 708 else (708.0-2*squareLB)/2.0

        for i in range(n_rows):
            for j in range(n_cols):
                it = self._scene.addRect(QtCore.QRectF(0,0,squareLB,squareLB),self.pen)
                it.setPos(p + j*squareLB, i*squareLB)

                self._squares.append(it)

        self.v.setScene(self._scene)
class Settings(QtGui.QMainWindow):
    def __init__(self, parent=None):
        super(Settings, self).__init__(parent)

        self.folder = QtGui.QPushButton("Folder", clicked=self.showSettings)
        central_widget = QtGui.QWidget()
        self.setCentralWidget(central_widget)
        vbox = QtGui.QVBoxLayout()
        vbox.addWidget(self.folder)
        self.scrollArea = QtGui.QScrollArea(widgetResizable=True)
        self.scrollArea.setBackgroundRole(QtGui.QPalette.Light)
        hlay = QtGui.QHBoxLayout(central_widget)
        hlay.addLayout(vbox)
        hlay.addWidget(self.scrollArea)
        self.setGeometry(200, 100, 300, 300)
    def showSettings(self):
        self.view = DemoApp()
        self.newwidget = QtGui.QWidget()
        hlay = QtGui.QHBoxLayout(self.newwidget)
        hlay.addWidget(self.view)
        self.scrollArea.setWidget(self.newwidget)
def main():
    app = QtGui.QApplication(sys.argv)
    ex = Settings()
    ex.show()
    sys.exit(app.exec_())
if __name__ == '__main__':
    main()


texts: How to set Persian charset in shapefile generated from GeoServer?

texts:
I have a .NET based web-gis application, and in that application different users update layers (layers information are stored in SQL server).

When I generate shapefile from that layer in GeoServer from WFS Service and then open it in ArcMap, ArcMap does not display the Persian values at all; instead it displays "????????".

In geoserver store i have used "windows1256".

I have tried different charset but I did not get the correct shapefile.


texts: Calculate the time spent by a user creating changeset

texts:
Is there a way to calculate the time spent by a user while creating a changeset?

The OSM Changeset file (i.e. map.osm) has timestamp, changeset, user
fields. Can the timestamp filed values be compared to find out
the total time spend by specific user?
Also, is there a way to do so using osmfilter?


texts: I already read this question Find EPS for DBSCAN but still I have some problem to find out what is the right value for 20 and 30 meters using WGS84 as reference for DBScan's eps. How can I find that?


texts: Calculating distance between a Point and a MultiPolygon PostGIS object

texts:
I'm working in a layer with MultiPolygons and I'm trying to calculate the shortest distance between a MultiPolygon and a point which is contained in that MultiPolygon, using Django and PostGIS.

At first I check which MultiPolygon contains that Point and I store it in a variable:

mpoly = MPoly.objects.filter(geom__contains=location)

Then, I'm trying to calculate the distance between the point and the mpoly object using the psycopg2 interface:

curs.execute("SELECT ST_Distance(ST_GeogFromText(%s), ST_GeogFromText('SRID=4326;Point(24.050, 35.311)')) FROM MPoly" % mpoly)

But it raises the following syntax error:

Traceback (most recent call last):
  File "<console>", line 1, in <module>
   ProgrammingError: syntax error at or near ";"
   LINE 1: SELECT ST_Distance(ST_GeomFromText(SRID=4326;MULTIPOLYGON ((23.567..

I think that the error comes from "mpoly"'s format but why ? If I type in a Python interpreter >> mpoly, it prints: MultiPolygon object at 0x7fc8d5130


texts: How to attach POIs in wikitude with google map places

texts: I want to create android application which uses AR and using camera it found the different stores around it and highlight it as display in this image. I am using wikitude sdk(javascript API) and android studio.I have created app like this image. I want attach all these POI markers with google map places.How to do that?


texts:
I have a few point shapefiles. I want to merge them, but first give each point a certain type according to shapefile. How do I add a String field, all with the same word?

(changing the field type from null to "water")

thanks.


texts: How to create one grid from specific set of points using QGIS?

texts:
I need to create one grid/fishnet from a set of points (coordinates). I tried to do it with the point to path tool following this three steps:

1-Vertical's rows: I numbered the points in the vertical of 1 to 8
       and created one vector with point to path tool.
2-Horizontal's rows: I numbered the points in the horizontal with odd
       numbers and created another vector with point to path tool.
3-Following, I merge this vectors.

However in the end result, I noticed that the lines created were parallel.

The lines pass in the points one more time. I tried separated the parallelal lines with vector --> data management tools --> split vector layer. I tried also vector--> geoprocessing tools --> dissolve. But both did not work. How can I solve this?


texts: ArcGIS Server 9.3 - Javascript API - Selecting Feature Problem

texts:
I have ArcGIS Server 9.3, and I'm trying to use the JS API. I successfully displayed the map, and I can locate and zoom to a feature based on ParcelID value. The problem is I can't highlight this parcel. I mean I tried to configure the symbol, but its not being displayed. Here's the code:

<script>

    dojo.require("esri.map");
    dojo.require("esri.layers.FeatureLayer");
    var map,featureLayer,query,sfs;

  function init() {
    map = new esri.Map("mapDiv", {
      sliderOrientation: "horizontal"
    });

    var imageParameters = new esri.layers.ImageParameters();
    imageParameters.format = "PNG24";  //set the image type to PNG24, note default is PNG8.

    //Takes a URL to a non cached map service.
    var dynamicMapServiceLayer = new esri.layers.ArcGISDynamicMapServiceLayer("http://*******/*****/rest/services/****/MapServer", {
      "opacity":1,
      "imageParameters":imageParameters
    });
    dynamicMapServiceLayer.setVisibleLayers([4,10,11,16]);
    map.addLayer(dynamicMapServiceLayer);

    //Locate Parcel--------
    sfs = new esri.symbol.SimpleFillSymbol(esri.symbol.SimpleFillSymbol.STYLE_SOLID,new esri.symbol.SimpleLineSymbol(esri.symbol.SimpleLineSymbol.STYLE_DASHDOT,new dojo.Color([255,0,0]), 2),new dojo.Color([255,255,0,0.25]));

    featureLayer  = new esri.layers.FeatureLayer("http://*****/*****/rest/services/*****/MapServer/4", {
            outFields: ["*"],
            mode: esri.layers.FeatureLayer.MODE_SELECTION
        });
    featureLayer.setSelectionSymbol(sfs);
    map.addLayers([featureLayer]);

    query = new esri.tasks.Query();
    query.where = "ParcelID_NS = '101010101'";

    var deferred = featureLayer.selectFeatures(query, esri.layers.FeatureLayer.SELECTION_NEW, function (selection) {
        var stateExtent = selection[0].geometry.getExtent().expand(3.0);
        map.setExtent(stateExtent);
    });
    //----------------------
  }

  dojo.ready(init);

</script>

Any help? Thanks.


texts: Geoserver configuration for large layers

texts: I'm having problem with geoserver configuration. Currently I'm running geoserver 2.4 on CentOS and Tomcat 6 with 24 cores CPU and 48 GB RAM. I have uploaded 4 layers to PostGIS (PostgreSQL 9.3/PostGIS running on the same machine) and when I want to sync those layers to geonode running geonode command manage.py updatelayers, after a few seconds geonode shows exception URLError: <urlopen error timed out>. While those errors are shown, all cpu cores are being used to 100%. I think that it's a problem with tomcat/geoserver configuration, could you give me some advice on how to set up tomcat/geoserver to host layers which contain more than 500,000 elements.


texts: How to USE ST_difference AND ST_intersection in case of multipolygons? (PostGIS)

texts:
I'm trying to use ST_difference and ST_intersection with two layers of multipolygons.

I want something like:

With the input on the left and the result on the right.

MY DATA:

t1.geom contains 1000 differents polygons.

t2.geom contains 1000 differents polygons.

THE QUERIES:

I try to create (in a naive way) two queries to perform these operations

st_intersection(t1.geom, t2.geom)
from t1,t2
where st_intersects(t1.geom,t2.geom);

where explain return:

"Nested Loop  (cost=0.15..49833.44 rows=149759 width=1796)"
"  ->  Seq Scan on affectation_com_1 t1  (cost=0.00..1447.85 rows=13085 width=923)"
"  ->  Index Scan using affectation_com_geom_idx on affectation_com t2  (cost=0.15..0.83 rows=1 width=873)"
"        Index Cond: (t1.geom && geom)"
"        Filter: _st_intersects(t1.geom, geom)"

AND

st_difference(t1.geom, t2.geom)
from t1,t2
where st_overlaps(t1.geom,t2.geom);

where explain return:

"Nested Loop  (cost=0.15..12768.08 rows=149759 width=1796)"
"  ->  Seq Scan on oca1032s_affectation_com_1 t1  (cost=0.00..1447.85 rows=13085 width=923)"
"  ->  Index Scan using oca1032s_affectation_com_geom_idx on oca1032s_affectation_com t2  (cost=0.15..0.83 rows=1 width=873)"
"        Index Cond: (t1.geom && geom)"
"        Filter: _st_overlaps(t1.geom, geom)"

Unfortunately it doesn't work, or takes forever to run.

I would like a query similar to the tools erase and intersect from ArcGIS.

QUESTION:

How to perform those operations efficiently ?


texts: ArcPy cannot open tables in Enterprise geodatabase

texts:
There are several scripts that create a geodatabase, modify the schema and insert data on the database. This scripts are designed as part of a solution. In a Nightly process this scripts are called, several times, in order to have several databases (Development/Testing process).
The scripts include the creation of SDE files for the connection based on User/Password.
When progress logs are checked,time to time, the scripts have failed:
cannot open 'DataOwner.PROJECT_STATUS'

but is not always the same table.
Running the set of scripts again in the morning they work fine.
Trying to find a possible cause, in the SQL server logs I found some entries that match the time when the scripts fail:
Date        8/11/2020 5:18:45 AM
Log     SQL Server (Current - 8/11/2020 6:00:00 AM)

Source      spid71

Message
The connection has been dropped because the principal that opened it subsequently assumed a new security context, and then tried to reset the connection under its impersonated security context. This scenario is not supported. See "Impersonation Overview" in Books Online.

Currently in the scripts, there is no impersonation.
After checking the configuration of the server, there is no other task/batch/script running at the same time that modifies/tries access to the database server instance. After some research with colleagues, the behavior has been seldom spotted on manual execution of the scripts.
What could be the possible causes for this behavior?


texts:
I'm trying to use an image overlay in KML / Google Earth as the basis for adding extra info and returning this to Google Earth.

In my utility app I can specify a rectangle to give me the north-south-east-west for the overlay and get the image, I then write this to a KML file as an overlay. (I can open this KML file in Google Earth and see that the overlay is nicely aligned).

I then import this KML file into my main app, which is a planar system. Thus I convert to UTM. The intention is to add info (i.e., draw on the image) and then export back from UTM to lat/long in another KML file. However, the 'rectangular' corners of the image in lat/long are not rectangular in UTM - so I need a transform in there some how.

How can this be achieved? Are the corners sufficient as GCPs? How do I calculate the necessary rotation/skew for the image? Is that even what I need to do?

(as you can tell, I'm no GIS expert!)


texts: Keyboard shortcut for saving and closing Query Builder window in QGIS 3

texts:
Pressing Enter while keyboard cursor are in the input box Provider Specific Filter Expression creates new line.
How to save and close Query Builder window in QGIS 3 without mouse?
This doesn't work:
Ctrl + Enter, Alt + Enter,
Esc, Alt + F4, Ctrl + S + Alt + F4  only closes the window, but does not save the results.


texts: ModisDownload error downloading only portion of hdf files

texts:
I am having problems with ModisDownload recently that I have not had when using it so far.
I dowloaded, reprojected and mosaiced modis images about two months ago successfully, despite some errors during downloading, which I assumed were due to the NASA servers, reported sluggish at times.

However, for the last week it is not working in a new way... I get no errors during downloading, but the .hdf files downloaded all are size 520 bytes, which is clearly wrong. MRT is then unable to open the .hdf files downloaded to reproject, leading to an error.

I have not changed anything to my configuration or my internet connection (from home, so no institution limitations on ftp access), and tried updating packages already.
Has anybody been experiencing the same error recently? Has something changed?

The code I have been using is below and is straightforward. I decomposed the historic download in loops because ModisDownload downloads first before reprojecting and mosaicing. I wanted avoid restarting from scratch every time there would be a downloading error.
I tried with a product name rather than the full http address, but I get the same error.

setwd(paste0("/media/olivier/olivier_ext/gedata_current/jde_coffee/MODIS/",state.nm,"/raw_data"))
FTP <- "http://e4ftl01.cr.usgs.gov/MOLT/MOD13Q1.005/"
years <- c(2015:2006)
start <- c("01.01","07.01")
end <- c("06.30","12.31")
for (i in 1:length(years)) {
  for (j in 1:length(start)) {
    ModisDownload(x=FTP,h=c(13,14),v=c(10,11),dates=c(paste0(years[i],".",start[j]),paste0(years[i],".",end[j])),
              mosaic=T,MRTpath='/home/olivier/MRT/bin', #'C:/modis_MRT/bin'
              UL=c(state_ex@xmin,state_ex@ymax),LR=c(state_ex@xmax,state_ex@ymin),
              bands_subset="1 1 1 1 1 1 1 0 0 0 1 0",
              proj=T,proj_type="GEO",proj_params="0 0 0 0 0 0 0 0 0 0 0 0 0 0 0",datum="WGS84",pixel_size=0.002259948,
              delete=T)
  }
}

I am right now switching to downloading by hand and then using reprojectHDF and mosaicHDF, but it is quite inconvenient for regular upadtes with new MODIS images.


texts:
I am new to the QGIS python console.
I am trying to load a raster layer from the console, set the crs and a specific style from a qml file.

layerPath = os.path.join(raster_dir,raster_name)
fileInfo = QtCore.QFileInfo(layerPath)
baseName = fileInfo.baseName()
layer = QgsRasterLayer(layerPath, baseName)

I would like to set up the projection info and then load it to QGIS

QgsMapLayerRegistry.instance().addMapLayer(layer)

I do have a grey raster on the windows. How can I set up the qml file to load ?

Thanks for your help


texts: Openlayers - Loop through control collection

texts: Is there a way to loop through the collection of map controls that have been added to an OpenLayers map?


texts:
I've imported a shapefile of the US and set a projection. I also added an Excel spreadsheet with income data from 1995-2010. I joined the income data with shapefile based on the state_name field.

What I would like to display is the percent increase or decrease over those ten years.  Is there a way to do this in Layer Properties | Symbology, or do I need to do the calculation in Excel before importing and joining?


texts: How to publish the Manage map server cache tile (Server) tool as an arcGis geoprocessing service?

texts:
I need to publish the  Manage map server cache tile on the arcgis server as a geoprocessing service. As a result I will be able to update the tiles of a tiled service through a service call.

I face 2 issues.
Fistly, when analyzing it tells me that I need to copy the datasource to the server (datasouce is the published tiled service used to RUN the tool once so as to be able to publish it).

Second it tell me that staging fails

Any ideas ?


texts:
I get the message "the plugin will be disabled because it crashed QGIS during last startup." Which plugin? None is identified. I'm a bit of a rookie so not sure where to look to solve this.
Running on a Mac with Mountain Lion, using QGIS Essen 2.14


texts:
I'm new to Carto and I'm trying yo upload a gdb file but I cannot see the option in the display formats.
How can I upload it?


texts: Converting vector to raster on-the-fly using GRASS?

texts:
EDIT: I later realized that, instead of creating multiple 256 x 256 tiles, I could just create a single image to overlay on the entire map. However, I'd still like to know if there are better ways to do the below.

I'm trying to create OpenStreetMap style "slippy" tiles (but mine are equirectangular, not Mercator projection) on the fly from a SHP file. These GRASS commands work:

# restrict the region to the slippy tile I want and set the rows/cols to
# create a 256 x 256 map tile

g.region n=41 s=40 w=-105 e=-104 rows=256 cols=256

# convert my existing vector map (which came from a SHP file) to a raster map
# Note: I reuse/overwrite 'temp' since I don't want to create dozens on
# unnecessary raster maps when I run commands like these several times

v.to.rast --overwrite input=ne_10m_time_zones output=temp use=cat

# convert the result raster map to a PNG file
r.out.gdal --overwrite input=temp output=somefile.png format=PNG

but I'm wondering if there's a better way. In particular, I'm wondering if this is even the right approach. Note that gdal_rasterize will do this from the command line but is much slower. In my approach, I only open the grass74 shell once and keep it open, which might be improving my time.

I know the standard method is to pre-create map tiles, but I'm interested in the on-the-fly methods as well.


texts: Running Create Extension postgis gives ERROR could not open extension control file?

texts:
I have installed PostgreSQL 9.x and PostGIS 1.5/2.0 several times and have never had this problem.

I just got a new CentOS 6.3 server up and running and and I have Postgres 9.3 functioning as expected. I have run

yum install postgis2_93

and I can see the files in

/usr/pgsql-9.3/share/contrib/

however, when I run

CREATE EXTENSION postgis;

I receive

ERROR:  could not open extension control file "/usr/pgsql-9.3/share/extension/postgis.control": No such file or directory

the tutorials (#1, #2) I have used don't show any steps between installing PostGIS and creating the extension.

What am I missing?


texts: Extracting raw data from WCS raster using QGIS python console?

texts:
I am trying to replicate the following manual workflow in the Python console of QGIS 3.4.7

Manual workflow:

1) Connect to WCS service.

Name: USGS National Map Elevation OGC-WCS

URL: https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer?request=GetCapabilities&service=WCS

2) Add the elevation data to the map canvas. I did this by dragging in the WCS layer "DEP3ElevationPrototype"

3) Zoomed to my area of interest

4) Right click the WCS layer > Export > Save As

5) Checked the box for Raw Data. Specified output file name and extent. Clicked OK.

Problem
The above workflow worked perfectly for me. However, I would like to have this same workflow in a python script. I couldn't find a processing tool for "Export > Save As", so I tried using GDAL Clip raster by mask layer

import processing
from qgis.core import *

### Raster layer from WCS
rlayername = 'elev'
uri = QgsDataSourceUri()
uri.setParam('url', 'https://elevation.nationalmap.gov/arcgis/services/3DEPElevation/ImageServer/WCSServer?request=GetCapabilities&service=WCS')
uri.setParam("identifier", rlayername)
rlayer = QgsRasterLayer(str(uri.encodedUri()), 'DEP3ElevationPrototype', 'wcs')

mask_filepath = 'N:\zTEMP\Ian\QGIS_Scripting\Black Diamon Study Area.shp'
mask_layer = QgsVectorLayer(mask_filepath, 'mask', 'ogr')

### Params for GDAL Clip Raster by Mask
parameters = {'INPUT': rlayer,
            'MASK': mask_layer,
            'NODATA': -9999,
            'ALPHA_BAND': False,
            'CROP_TO_CUTLINE': True,
            'KEEP_RESOLUTION': True,
            'OPTIONS': None,
            'DATA_TYPE': 0,
            'OUTPUT': 'N:\zTEMP\Ian\QGIS_Scripting\DEP3_Extract_Clip.tif'}

### Clip
print ("running algorithm")
processing.runAndLoadResults("gdal:cliprasterbymasklayer",parameters)
print("done")

I assume the problem is that Clip Raster by Mask cannot take a WCS raster, though I could be mistaken. Are there any ways to doing something like the manual workflow in the Python console?

Note that in this case I don't really care about the distinction clipping by extent vs clipping by mask. In the end I just want to get a TIF extract from a WCS raster.


texts: What are the benefits of normalizing longitude to [-180, 180)?

texts:
In the Google Maps API v2 for Android, LatLng normalizes longitude values to the range -180 degrees inclusive to +180 degrees exclusive.  Consequently, attempting to represent the entire world as a LatLngBounds produces a longitude range of -180 to -180, which has zero width and contains no points:

// southwest, northeast
final LatLngBounds world = new LatLngBounds(
    new LatLng(-90, -180),
    new LatLng(90, 180));
assertTrue(world.contains(new LatLng(0, 0))); // test fails

So obviously, normalizing longitude to [-180, 180) has a cost in terms of how LatLngBounds works.  We can only speculate as to why the API designers might have made this decision.  But what are the potential benefits of normalizing longitude this way?  Are there any common operations that are simplified?


texts: Where can I find city area data for the US?

texts: Where can I find (on the site of a US government agency, I assume) raw data for the sizes of cities that I can programmatically process? Data indexed by FIPS place codes (i.e. those found here) would be optimal.


texts: How can I enter metadata into a qgiscloud project

texts:
I published a project into qgiscloud: https://qgiscloud.com/tragewegen/IDEAproject_DonkeyMobileBlackCarbon.

I want metadata to be published in the metadatatab next to the legendtab.  Can somebody tell me how to do that?


texts: Projecting 2D (poly)lines and polygons on point cloud or (tri)mesh (TIN)?

texts:
I'm searching a way to project 2D (break)lines and polygons (they are in an Esri shapefile but I can translate them to something different if needed, like geojson for example) on a point cloud (seems a priori harder because of the definition of what is a projection on a point cloud) or a triangular mesh (*.ply format, or eventually other formats if you have better proposals) and thus having new points along the triangles edges with Z attributes.

These 2D lines need to have a new Z attribute (in geom or attribute table, I don't know if 3D shapes does exist) with altitude after the operation.
And the mesh should also have these new points.

I'm using open sources softwares under Ubuntu like PDAL, PCL, CloudCompare, QGIS, SAGA-GIS and libraries (GDAL/OGR).
GIS are not always nice with huge amount of 3D data.
Command line interface-tools are welcome.

About loading point file formats:

first: with *.ply files : QGIS can't open them. So I was trying to convert these in a vector format (*.shp) using for example CloudCompare (it's for the moment the only software I found that can do that). But it only keeps the vertices (I can't find a way to import a 3D triangular mesh, for example by importing the vertices in a shape and the triangles in an other, which would be a nice way to do that), not the edges of triangles from the original mesh *.ply, so I will have to re-run a triangulation algorithm (see below) within QGIS.
Then: with *.las files; using lastools' las2shp or the same export to *.shp process as before will also convert my points of interests to a vector format (shp). But I will then also need to run a triangulation algorithm. Prior to that, I've tried to mesh the *.las and then exporting the triangles vertices as explained for the *.ply format.

Then, with the rendered points in vector format:

first: it's really slow (display and math on them), it take about minutes only to refresh the zoom level (I fear a little about math or SQL operations on them...)
Then: I'm stuck with the triangulation process. As some vertices are "under" others, they shouldn't be all used, i.e. all triangles faces should be oriented upwards because, AFAIK, GIS softwares can't handle true 3D; they uses 2D Delaunay. Geometric intersections are also only 2D in GIS softwares. So, if I intersect 2D lines with the triangles edges, I may have "errors" with edges that are below other and that shouldn't exist topologically.
This is because two points that seem OK to be linked in 2D by the (2D)Delaunay triangulation should not necessarily be linked (some surely were not, for example in the previous 3D mesh generation process, outside GIS, e.g. in CC or meshlab) as they are really far each others in 3D (but not in 2D)!

After that I can't see how to get the Z value on an intersection between a triangle edge and my 2D line of interest which is really 2D and that I want to project on the triangles. Projecting, here, is of course the same as extracting the Z values from the triangles' edges.


texts: 404 config.json error for server deployment of ArcGIS Web AppBuilder Application?

texts:
I am receving  a 404 error that my config.json cannot be found, I have added the JSON mime type to my IIS configuration at the site and server levels and have restarted IIS.  This error appears as so in my console  with other errors which I think are caused by the initial 404 of the config.json.  What other steps can be taken to remedy this?  Additionally when I hit the index.html I receive  an "Unexpected token < in JSON at position 0" error

GET http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/css/fonts/c4352a95-7a41-48c1-83ce-d8ffd2a3b118.woff
main.js?wab_dv=2.1:1521 jimu.js init...
init.js:140 GET http://gisdev01/LifeSupportApplication/LifeSupportApplication/config.json?wab_dv=2.1 404 (Not Found)h @ init.js:140_tryLoadConfig @ main.js?wab_dv=2.1:236loadConfig @ main.js?wab_dv=2.1:227loadConfig @ main.js?wab_dv=2.1:11f.initApp @ main.js?wab_dv=2.1:1522(anonymous function) @ init.js:200ia @ init.js:28(anonymous function) @ init.js:28ja @ init.js:28ga @ init.js:28b @ init.js:29f @ init.js:151d @ init.js:152
init.js:113 TypeError: Cannot read property 'portalUrl' of undefined(…) "TypeError: Cannot read property 'portalUrl' of undefined
    at _addDefaultPortalUrl (http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/main.js?wab_dv=2.1:32:492)
    at _addDefaultValues (http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/main.js?wab_dv=2.1:31:416)
    at .<anonymous> (http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/main.js?wab_dv=2.1:11:456)
    at http://js.arcgis.com/3.17/init.js:63:209
    at h (http://js.arcgis.com/3.17/init.js:107:279)
    at k (http://js.arcgis.com/3.17/init.js:107:205)
    at g.resolve (http://js.arcgis.com/3.17/init.js:109:280)
    at a (http://js.arcgis.com/3.17/init.js:108:148)
    at h (http://js.arcgis.com/3.17/init.js:107:419)
    at k (http://js.arcgis.com/3.17/init.js:107:205)
    ----------------------------------------
    rejected at a (http://js.arcgis.com/3.17/init.js:108:174)
    at h (http://js.arcgis.com/3.17/init.js:107:438)
    at k (http://js.arcgis.com/3.17/init.js:107:205)
    at g.resolve (http://js.arcgis.com/3.17/init.js:109:280)
    at a (http://js.arcgis.com/3.17/init.js:108:148)
    at h (http://js.arcgis.com/3.17/init.js:107:419)
    at k (http://js.arcgis.com/3.17/init.js:107:205)
    at g.w.reject (http://js.arcgis.com/3.17/init.js:109:488)
    at a (http://js.arcgis.com/3.17/init.js:108:174)
    at h (http://js.arcgis.com/3.17/init.js:107:452)
    ----------------------------------------
Error
    at then.c.then (http://js.arcgis.com/3.17/init.js:110:97)
    at loadConfig (http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/main.js?wab_dv=2.1:11:361)
    at Object.f.initApp (http://gisdev01/LifeSupportApplication/LifeSupportApplication/jimu.js/main.js?wab_dv=2.1:1522:55)
    at http://gisdev01/LifeSupportApplication/LifeSupportApplication/init.js:200:24
    at ia (http://js.arcgis.com/3.17/init.js:28:3)
    at http://js.arcgis.com/3.17/init.js:28:256
    at ja (http://js.arcgis.com/3.17/init.js:28:125)
    at ga (http://js.arcgis.com/3.17/init.js:28:198)
    at b (http://js.arcgis.com/3.17/init.js:29:177)
    at f (http://js.arcgis.com/3.17/init.js:151:195)"m @ init.js:113(anonymous function) @ init.js:114filter @ init.js:71h @ init.js:114


texts: Precision of x decimal points in meters?

texts:
Is there a handy reference that corelates the number of decimal points in a coordinate to the precision in meters?

I had discussions like "4 decimal places will place you within a 500 meter radius" and "6 decimal places will place you in a house". I want a reference that can clear up such facts and requirements.


texts:
We're in the process of transforming a lot of 9.1 ArcServer/ArcObjects etc applications to ArcGIS Server 10.

However, we're unsure as to licensing on 10.

As I understand it, I cannot, according to license restrictions, deploy a toolbox that has been written using ArcObjects 10 .NET SDK on the ArcGIS 10 Server for Java platform, for licensing issues.  We can deploy a .NET application to an ArcGIS Java Server 10 instance, if we bundle up all of the application object references in the installer, but this goes against the license agreements (you're not allowed to bundle them up in the installer apparently), so these have to be deployed on an ArcGIS 10 .NET Server instance.  Is that pretty much the long and the short of it?

I am sure we can't have two different flavour instances of ArcGIS 10 Server installed on the same machine (one Java, one .NET) so they need to be on separate boxes (or VM's)?

Thanks in advance for any words of calm...


texts:
I am supporting a .NET code that is running on ArcGIS Server 10.1, and communicating with ArcSDE Application Server 9.3

There is part of the code that is doing long loop.
And when it comes back to communicate to the database (to ArcSDE), it is raising this exception
ORA-03113: end-of-file on communication channel

It seems a timeout error.
Anybody has any idea how to fix this?


texts: ArcGIS Desktop: Dynamic label of layer attribute (Sum of line length)

texts:
I'm trying to put a dynamic text label beside a legend entry for a road layer that has the length of the road dataset in ArcMap Desktop (10.6).
I usually work in QGIS so I'm not familiar with the proper way to do this in ArcMap. In QGIS you could create a label and use this code:
[%to_int(aggregate('my_layer','sum',"FEATURE_LENGTH")/1000)%]
However, I'm not sure what the equivalent would be in ArcMap and can't seem to dig up an answer with Google. Does anyone know how to do this?


texts: Making Hibernate Spatial DWITHIN compatible with both H2GIS and PostGIS

texts:
I have following:
@Query("SELECT f FROM Feature f WHERE DWITHIN(f.geometry, :centre, :range, true) = TRUE")

DWITHIN is hibernate spatial function explained here.
boolean dwithin(Geometry, Geometry, double) Returns true if the geometries are within the specified
distance of one another

In my query above I pass 4 parameters. Internally hibernate spatial maps it to postgis ST_DWithin which has 2 variants. If I don't pass fourth parameter, first variant is considered  which takes 3rd parameter as SRID degree range and not distance. If I pass fourth parameter, third parameter is considered as distance in meters.
But for H2GIS(OrbisGIS), there is only a single variant which takes 3 parameters, and it always takes 3rd parameter as distance.
For above query to work with both, I would need to get rid of fourth parameter, in which case PostGIS version would start treating 3rd parameter as SRID degree range while I want it to be treated as distance in meters.
This is my understanding of the issue.
Am I right here? If yes, how can I make it work for both, accepting 3rd parameter as distance in meters.


texts: Find highest point of 3D line with PostGIS

texts:
I need to find the highest point of a 3D line. To do so I use ST_ZMax to obtain the Z-coordinate of the highest point:

SELECT ST_ZMax(geom) FROM 3dlines;

How can I find the X- and Y-coordinates of this point?