texts: Define categorized renderer with symbols based on several categories

texts:
In pyQGIS I have added a vector layer based on a CSV file. The CSV file contains several fields: x , y, Capacity, Category.

I can render the map with simple points at the given (x,y) coordinates.

What I would like to do give each point a circular symbol with a diameter based on the Capacity value and a color based on the Category value. This is my current code for changing the color:

project = QgsProject.instance()
project.read("filename.qgz")
uri='file:///C:/Till//csv_file.csv?delimiter=,&yField=y&xField=x'

new_layer= QgsVectorLayer(uri, "new", "delimitedtext")
project.addMapLayer(new_layer)

fni = new_layer.fields().indexFromName('Category')
unique_cats = new_layer.uniqueValues(fni)
categories=[]

for unique_value in unique_cats:
    symbol = QgsSymbol.defaultSymbol(new_layer.geometryType())

    layer_style = {}
    layer_style['color'] = '100, 100, 100' # just for testing
    layer_style['outline'] = '#000000'
    symbol_layer = QgsSimpleFillSymbolLayer.create(layer_style)

    if symbol_layer is not None:
        symbol.changeSymbolLayer(0, symbol_layer)

    # create renderer object
    category = QgsRendererCategory(unique_value, symbol, str(unique_value))
    # entry for the list of category items
    categories.append(category)

renderer = QgsCategorizedSymbolRenderer('EdgeRef', categories)

new_layer.setRenderer(renderer)

However when I render this layer, it just shows empty. What is going wrong? And how can I add the diameter feature to the geometry?


texts: Add a local GeoJSON file as a layer with GeoJsonLayer

texts:
I turned a state boundary feature class into a geojson file called region2_boundary.json:

 {
  "type" : "FeatureCollection",
  "crs" : {
    "type" : "name",
    "properties" : {
      "name" : "ESRI:102008"
    }
  },
  "features" : [
    {
      "type" : "Feature",
      "id" : 1,
      "geometry" : {
        "type" : "Polygon",
        "coordinates" : [
          [
            [
              1708004.4871000014,
              329793.41640000045
            ],
           ....
  }

I have placed it into my Scripts folder in my ASP.Net app at the same level as my app.js. I am trying to bring it in as a layer in my app.js code, but there seems to be a problem. I keep getting the error:

  Uncaught ReferenceError: GeoJsonLayer is not defined

I am thinking that I am just using a bad filepath/url in my boundaries variable. Any suggestions on how to get this simple json file added as a layer or how to resolve this error?

require([
  "esri/map",
  "esri/basemaps",
  "esri/dijit/BasemapToggle",
  "esri/layers/ArcGISDynamicMapServiceLayer",
  "esri/tasks/ImageServiceIdentifyTask",
  "esri/tasks/ImageServiceIdentifyParameters",
  "esri/layers/WebTiledLayer",
  "esri/dijit/Search",
  "esri/config",
  "esri/request",
  "esri/tasks/IdentifyTask",
  "esri/tasks/IdentifyParameters",
  "esri/tasks/query",
  "esri/tasks/QueryTask",
  "dojo/json",
  "dojo/promise/all",
  "dojo/domReady!"], function(
    Map,
    esriBasemaps,
    BasemapToggle,
    ArcGISDynamicMapServiceLayer,
    ImageServiceIdentifyTask,
    ImageServiceIdentifyParameters,
    WebTiledLayer,
    Search,
    esriConfig,
    esriRequest,
    IdentifyTask,
    IdentifyParameters,
    Query,
    QueryTask,
    JSON,
    all) {
    ...
    var boundaries = "./region2_boundary.json";
    var geoJsonLayer1 = new GeoJsonLayer({
        data: boundaries
    });
    ...
    map = new esri.Map("map", {
        //basemap: "streets",
        basemap: "topo",
        center: [-66.664513,18.200178], //PR
        zoom: 9,
        layers: [geoJsonLayer1]
    });

EDIT I have added my dependencies above.


texts: write multiple statistics analysis into the same output file with arcpy

texts:
How can I add multiple stats from the arcpy.Statistics_analysis to a single output file. Can I create a text file instead of a table ?

Using one of the examples in help.arcgis, if I add other stat analysis, what is the best way to work with the output?

import arcpy
from arcpy import env
env.workspace = "C:/data/Habitat_Analysis.gdb"
arcpy.Statistics_analysis("futrds", "C:/output/output.gdb/stats1", [["Shape_Length", "SUM"]], "NM")
arcpy.Statistics_analysis("futrds", "C:/output/output.gdb/stats2", [["Shape_Length", "MEAN"]], "NM")
arcpy.Statistics_analysis("futrds", "C:/output/output.gdb/stats3", [["Shape_Length", "STD"]], "NM")


texts: Converting distance between nautical miles and degrees using ArcObjects?

texts:
We have a certain functionality developed that:

Will check for a tolerance value (this input is given by the user as NM);
Buffer a polygon by that value (all polygons are in WGS-84);
Iterate through selection (buffer x points) and creates a record for each intersected point.

The tester at my corporation is saying the tool is giving wrong results. Here is his workflow, for testing:

Creates a polygon;
Define a nearby point;
Measure the distance between the polygon and point (using measure tool);
Configure tolerance to be larger or equal distance;
RUn the tool;

I have a problem in the middle, is that our user needs to store these tolerance values as nautical miles, and I need to convert them to degrees for them to be used in our buffer (I'm using IUnitConverter). All buffers are constructed using ITopologicalOperator.

My guess is the conversion is being done purely in a numerical way and then, buffering by less than it should.

I've read tons of articles speaking about projection engine.dll, but that is for calculating distances, not buffers (I mean, if I need to implement a custom buffering solution, that will take into account each vertex).

Can you guys recommend an approach for this? How can I convert this measure and still hit these edge cases?


texts:
I need to visualize in QGIS the density of photographs taken in seven municipalities of the central Italy: Perugia, Magione, Passignano, Umbertide, Deruta, Corciano, Torgiano.

I'm wondering if there is a method to export the Panoramio georeferenced features (e.g. the places where photos were taken) in a format that can be opened in QGIS (csv, shp or others), or if there is a method to import the Panoramio georeferenced features as OpenLayers in QGIS.

Can anybody help me?


texts: How to write a expression in field calculator for finding $X and $Y?

texts: I'm trying to calculate the lat/long of points using the Q GIS field calculator.Nathan W suggested to First, save the layer as WGS84 and import it again. Then in the field calculator you can use $x or $y as a variable to get the x and y. but the problem is this: i dont know how to write a expression in field calculator for finding $X and $Y?


texts: How can i simply merge csv files together with fme

texts: I have some csv files which have varying contents. I want to simply disregard the schema and write the lines of text into one file. Which transformers or technique could I use to achieve this?


texts:
During the tutorial : http://www.qgistutorials.com/en/docs/performing_spatial_queries.html

I have a error message like below.

[

What is the providers mean?

I really try to figure it out but nothing I can find.

My qgis version is QGIS 2.18.19


texts:
I'm working on an application with a lot of GPS sensors, for each vehicle and person in a scene.

I want to show a vehicle and his occupants in the same marker. But I can't get the vehicle occupants list from the system, I only have GPS positions and object type (firetruck, policeman, car driver, etc…).

I'm looking for an existing model or library doing this kind of merge.

The simple approach can be to merge closest sensors from a vehicle into it, but sometimes you have a lot of people outside of a vehicle.

I can use the historical data, if some sensors are moving at the same speed and are near since a long time, I can consider that they are in the same vehicle. But it's a homemade system and it's maybe not very good.

Is there a good way?


texts:
i have problem with my json of delete feature form Applyedit Jason on the arcgis server Json. Here my code.

delMaintData:  function(){
console.log('Bri');
            var maintTable = new FeatureLayer(this.dataSelect.get("value"));
            console.log('layer',maintTable);
            var deletesrecord = this.OBJECTID.get('value');
            console.log(deletesrecord);
            maintTable.applyEdits( null, null, [deletesrecord], function(deleteResults){
                //readTable(dijit.byId('assetId').get('value'));
                //dijit.byId('dialogAsset').hide();
                console.log(deleteResults);
            });

    },

i have [deletesrecord] to define the object id and i get the number          from console.log(deletesrecord); , and after the code in maintTable.applyEdits , there are log error

Uncaught TypeError: Cannot read property 'undefined' of undefined


texts: Bind Vector Geometry to ViewModel or How to keep track of modified points?

texts:
I'm trying to create a viewModel for my application and the data is typically populated by a call to a database. I can retrieve the coordinates and convert this to geometry which can be displayed on a map but I would like to bind these two together so that if a feature is modified (say points are added or removed) the information is relayed back to the viewModel for updating to the database.

The main difficulty I've encountered is that I can't set properties on the individual vertices of a linestring so trying to get the ID of an individual vertex and then finding the corresponding id on the VM isn't an option.

Alternatively I can try to get the index of the item in the two respective arrays and use an array method but that just seems entirely hit-and-miss. I've been considering utilising KnockoutJS in order to implement some rudimentary binding but I'm wondering if there isn't something baked in to OL3 by way of the ol.Observable class.

More importantly, if I wished to show rich data for each point such as elevation, name, etc., how could this possibly be attached to the geometry point itself? I've encountered a jsFiddle where each vertex of a linestring is actually a clone of a Point object and moving one will be reflected in the other but it seems like it will take an enormous amount of effort for me to re-engineer how my drawing interactions work. Can anyone enlighten me as to how this might be achieved?


texts: Getting a plugin path using Python in QGIS

texts:
How to get a plugin path using Python?  I tried QgsApplication.pluginPath()

I got u'/Applications/QGIS.app/Contents/MacOS/../PlugIns/qgis'.

But I actually store the plugin in ~/.qgis2/python/plugins/myplugin.


texts: QGIS processing.run() standalone script does not write .prj and .cpg output files of shapefile?

texts:
I've created a python standalone script to clip two shapefiles with qgis. I've installed qgis3 with OSgeo4W. The logic is as follows:

import qgis
from qgis.core import *
from qgis.core import QgsApplication
from qgis.core import QgsProject
from qgis.core import QgsProcessingFeedback
from qgis.core import QgsVectorLayer
from qgis.analysis import QgsNativeAlgorithms
import processing
from processing.core.Processing import Processing

path_qgis_installation = "C:/OSGeo4W64/apps/qgis"
QgsApplication.setPrefixPath(path_qgis_installation, True)

qgs = QgsApplication([], False)
qgs.initQgis()

Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

input_file = QgsVectorLayer('.../file_name.shp')
overlay_file = QgsVectorLayer('.../file_name.shp')
output_fil = '.../file_name.shp'

processing.run(
    "native:clip",
    {
        'INPUT':input_file,
        'OVERLAY':overlay_file,
        'OUTPUT':output_fil
    }
)

qgs.exitQgis()

My problem is that when I run the command (only the processing.run()) in the python console in qgis or in the osgeo shell, everything works.

When I however run it with my python console (bat script), I only get the following (corrupted) files written to disc: .dbf, .shp, .shp. Missing are .prj and .shx for a complete shapefile.

My suspition is that is has something to do with the paths. At one point I executed C:/OSGeo4W64/bin/o4w_env.bat from a location outside of C:/OSGeo464) which then even procued the same problem whilst running on the osgeo shell.

Any suggestion what I am missing to load?


texts:
I'm having issues with some scripts that import an arc editor licence to run (similar to the example given on this page http://resources.arcgis.com/en/help/main/10.1/index.html#//002z0000000z000000), yet do not return the licence to the licence manager once the script has been successfully completed.

Does anyone know of similar functionality to "arcpy.CheckInExtension, arcpy.CheckOutExtension", which could be applied to Product Modules (arcinfo, arceditor, arcview and so on)?


texts:
Is there a way to count the number of data driven pages a polygon is on in an ArcGIS mapbook?  I need to identify parcels that expand across multiple sections (which is my index for creating data driven pages), and I do not have any keys that would allow me to tie parcels to sections

I was able to write something up, the issue is I keep getting runtime errors because I'm looping a 50,000 item list through a 5,000 item list.  Any algorithm suggestions?

def get_extent(section_class, parcel_class, sc_fields, pc_fields):
    sc_cursor = arcpy.da.SearchCursor(section_class, sc_fields)
    pc_cursor = arcpy.da.SearchCursor(parcel_class, pc_fields)
    poly_count = ['']
    for sc_row in sc_cursor:
        sc_extent = sc_row[0].extent
        for pc_row in pc_cursor:
            print("Working on: ", pc_row[1])
            pc_extent = pc_row[0].extent
            if fc_extent.crosses(pc_extent):
                print("Parcel {0}cross map {1}".format(pc_row[1], fc_row[1]))
                poly_count.append(pc_row[1])
    print("All done, here's the list: ", poly_count)
    return poly_count


texts: How to convert GML to features in ArcObjects?

texts:
How to convert GML to features and features to GML i ArcObjects?
We are building an application where we want to communicate with a service with GML but natively in our stand-alone application we want to use features and shapefiles.

Update: I have an ArcGIS Server so I could use that a bridge so I could use thatas a bridge beween the service and the client.


texts:
I've been searching for quite a while now. Is there a way to merge a hillshade image with an aerial image with gdal?
I'm just trying to create a hillshaded aerial image.


texts: PostgreSQL PostGIS ST_Point_Inside_Circle() radius in meters

texts:
I'd like to begin by saying that I have very little knowledge and understanding when it comes to GIS (despite extensive experience with Google Maps API). PostGIS appears to be a whole new ball game on its own.

I've recently installed PostGIS (v2.1.5) with my PostgreSQL (v9.4) instance in order to begin storing and working with GIS data. I'm currently using Google Maps API v3's DrawingManager to draw polygons, circles, rectangles and markers which need to be stored in the database, and later retrieved to display the polygons again on Google Maps; These shapes will also be used for various aggregate queries.

I currently have two tables. The one is called "locations" which will store the polygons, circles and rectangles created with Google Maps' DrawingManager. The second, called "rawdata" will contain a set of records with a POINT (lat, long) inside of it.

First and foremost, I do not entirely understand the concept of projections and how they work, but, based on what I've read online is that calculations and accuracy may vary, significantly. I understand that Google Maps works with the Web Mercator projection which I believe is officially known as the EPSG:3857 projection (or CRS), while the preferred method for storing and working with geometry() shapes in PostGIS is in the EPSG:4326. As such and in reference to my create table statements below, I have chosen to store my geometry data in EPSG:4326 in the tables.

An easier method would to be store my geometry in the newer geographic format, but I do not wish to impose those performance and lack of functionality drawbacks as I will be working with a significant amount of data (hundreds of thousands, if not millions of records).

I am trying to test the ST_Point_Inside_Circle() function but this appears to be inaccurate and I can only assume so because the function expects degrees as the radius parameter, yet Google Maps supplies me with a radius in meters. In reference to http://resources.arcgis.com/en/help/main/10.1/index.html#//009t0000004v000000, I have tried converting meters to degrees by multiplying meters by 0.00001. However, the greater the radius in meters, the more inaccurate it appears to be. What can I do to achieve an accurate result?

SELECT ST_Point_Inside_Circle(ST_SetSRID(ST_MakePoint(28.239418677985668, -25.75443387513752), 4326), 28.241967111825943, -25.759674774140656, (583 * 0.00001));

Google Maps reports the distance in meters between the two points (center of circle and point) as 637 meters, yet ST_Point_Inside_Circle returns true when trying to convert 583 to degrees. This leads me to the assumption that one cannot convert meters to degrees using 0.00001.

What is the correct way to determine the radius with meters?


texts:
I'm trying to install PyWPS on my computer and when I run

 sudo python setup.py install

I'm receiving :
"extensions/gdal_wrap.cpp:2853:22: fatal error: cpl_port.h: No such file or directory
compilation terminated.
error: Setup script exited with error: command 'x86_64-linux-gnu-gcc' failed with exit status 1"

I've Gdal 2.1.2 installed on my computer and I can see the cpl_port header in /usr/include/gdal/.

Is there something I'm missing ?


texts: Cross information in ArcGIS for Desktop?

texts:
I'm a total amateur at ArcGIS for Desktop, this is the first time I use it for real.
The thing is I have several shapefiles with information about:

Schools in my city (elementary, middle and high school)
City districts

The schools' information has Points that represents the unit (the place itself) and also a Buffer (influence's radio of every unit). ONE LAYER BUT SEPARATE SHAPES (One for point one for buffers).

So, I need to cross this information so I can generate data that tells me the districts that are best equipped with each kind of educational level. I mean, which districts have more coverage of elementary schools, more coverage of middle high schools and more coverage of high schools. Even if they don't have any school within, but just some coverage from any buffer.

I don't know if I explained myself very well, but I'll be so thankful for any advice.
I think I must merge the schools and then join this shape with the districts... (From what I've learned in my two days ArcGIS experience)


texts: QGIS georeferencer: Thin Plate Spline (TPS) residual

texts:
1) Is anyone able to explain in simple manner what does mean error exactly mean? I use metric values, so suppose 3.64109 is expression of meters and "e" is mathematical constant:
http://en.wikipedia.org/wiki/E_%28mathematical_constant%29
but what's next?

2) Why does image become more brighter after TPS transformation?

I work on archival aerial imagery and I'm testing all of the transformation types and re-sampling methods to find the best one. As you can anticipate, I'm not mathematician and I'd like to  understand what's result of my work in fact.


texts: Create a report wizard fails in ArcGIS 10 (win7 x64) - but works in ArcGIS 10 (winxp 32bit)

texts:
The report creation wizard will not launch for me when I try accessing it from both the View menu and from an attribute table.  I am able to open a saved report template though.  This seems like it is a Windows 7 issue since the wizard does launch under Windows XP.  Or it could be a 64bit vs 32bit issue.

This has been reported in the ArcGIS blog Working with Reports in ArcMap, Part I-Getting Started.

I actually ended up doing the task using Mail Merge in Word 2007.  It was so easy to format it together from an exported attribute table - I'd definitely do this again, even if I could get the report wizard to work.


texts: From uploaded zipped shapefile to Geopandas DataFrame in Django Application

texts:
I'm trying to get to a point where I can quickly filter thousands of points in a shapefile. My Django application asks for a zipped shapefile to upload, where the zipped file contains at least the .shp, .shx, and .dbf files. Once in my Django view, the zip file is as follows:
request.FILES['file'] > <InMemoryUploadedFile: test.zip (application/x-zip-compressed)>
type(request.FILES['file']) > <class 'django.core.files.uploadedfile.InMemoryUploadedFile'>
request.FILES['file'].file > <_io.BytesIO object at 0x0000028E29F8FE00>
Assuming Geopandas is the best option for efficient filtering/masking (if I'm wrong, I'm definitely open to suggestions), I'm not sure how to go from current state to a Geopandas DataFrame. When I try to use the read_file() method
import geopandas as gpd
gpd.read_file(request.FILES['file'].file)

I get the following error:

fiona.errors.DriverError: no driver

The geopandas.read_file() docs state:

Either the absolute or relative path to the file or URL to be opened,
or any object with a read() method (such as an open file or StringIO)

I'm not sure how to get what I have into an appropriate format for the read_file() method.
Note: The masking and filtering I'm looking to perform are on attribute data and not the geometry.


texts: Adaptive Buffer Distance Expression in QGIS

texts:
I have a shapefile of line features that I need to add a buffer to. The width of the buffer will depend on the length of each feature (25% of the feature length). I know I need to create an expression for the buffer distance, but have not been able to come up with the proper expression.  Any idea how to change the expression to come up with the right result?

I have tried the suggestion shared here but the buffers still end up the same width for all of the features.


texts: Comparing and find out which feature has highest shape area in QGIS?

texts:
I need to find out which polygon is having more area by filtering unique ID. the problem is I have 1000 of features which I need to find out highest occupied area for every unique value...

since I have 1000 of features I need to do in an automated way instead of finding out manually,  please help me to solve with some Python script or by using field calculator queries.


texts: Insert point into point column in postgres/postgis

texts:
I'm trying to insert point data in my database, but it tells me the data types aren't correct (and I think they are). Here's my table:

 Column      |  Type
 ------------+----------
 coordinates | point

And my query:

INSERT INTO points(coordinates) VALUES (ST_GeomFromText('POINT(10.809003 54.097834)',4326));
ERROR:  column "coordinates" is of type point but expression is of type geometry
LINE 1: INSERT INTO points(coordinates) VALUES (ST_GeomFromText('POI...
                                                ^
HINT:  You will need to rewrite or cast the expression.

I don't quite get what is wrong (using the same as a select works). I get the same with LINESTRING and POINT, and my syntax seems to be the same as How to insert a point into postgis?

I partially made it work with help from @amball by dropping the table and re-making it. The insert worked, but not I have this:

gisdb=# INSERT INTO points(coordinates) VALUES (ST_GeomFromText('POINT(10.809003 54.097834)',4326));
INSERT 0 1
gisdb=# select * from points;
 id |                    coordinates
----+----------------------------------------------------
  1 | 0101000020E610000056BB26A4359E2540EB3713D3850C4B40
(1 row)

I don't know what happened with the coordinates (they are supposed to be 10.809003 54.097834 or something similar.

Here's the table info, which looks better, by the way:

   Column    |         Type         |                      Modifiers                      | Storage | Stats target | Description
-------------+----------------------+-----------------------------------------------------+---------+--------------+-------------
 id          | integer              | not null default nextval('points_id_seq'::regclass) | plain   |              |
 coordinates | geometry(Point,4326) |                                                     | main    |              |


texts: Pansharpen landsat mosaic in Google Earth Engine

texts: Does anyone have a good way to pansharpen an entire mosaic of landsat 8 images? I would use the Panchromatic band (band 8). I have found tutorials on how to pansharpen a single image, but not for a mosaic (or, essentially, an entire ImageCollection).Preferably, I would carry this out in Google Earth Engine.


texts: QGIS on Samsung GT-N5110 - invalid Qt version

texts: I want to install QGIS Experimental on my Samsung GT-N5110 tablet. When opening QGIS cannot open additional libraries and there is an error message: "Invalid Qt version" and QGIS crashes. I use Android Version 4.4.2. Are there any ideas to solve this problem?


texts: Converting Tiles to GeoPackage with Python

texts:
I am trying to convert tiles (as PNG) into GeoPackage. In order to do so, I would like to use the geopackage-python tool. However, when executing the script I run into the following error:

Layer 14 is the folder with the tiles/images.

Do you have any idea on what the problem exactly is and how to solve it?


texts: How to dynamically update a KML layer in OpenLayers 2.10?

texts:
Arguably, this is a easy question for you guys GIS expert. Here is the question, how to outsource a dynamic KML file using JavaScript in OpenLayers. I want to refresh every 7 seconds.

Below are snippet of the given codes:

<script type='text/javascript'>

var map;

OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3;

function init() {
    map = new OpenLayers.Map({
        allOverlays: true,
        controls: [
            new OpenLayers.Control.Navigation(),
            new OpenLayers.Control.PanZoom(),
        ]
    });

    var osm_layer = new OpenLayers.Layer.OSM();
    map.addLayer(osm_layer);

    var MyKML = "http://www.xyz.com/tld/mypath/myfile.kml"

    map.addLayer(new OpenLayers.Layer.GML("KML", MyKML,
           {
            format: OpenLayers.Format.KML,
            formatOptions: {
              extractStyles: true,
              extractAttributes: true,
              maxDepth: 1
            }
           }));

    // Adding controls to map
    map.addControl(new OpenLayers.Control.ScaleLine());
    map.render('map_element');

    map.setCenter(
            new OpenLayers.LonLat(101.587008,3.021612).transform(
                new OpenLayers.Projection("EPSG:4326"),
                map.getProjectionObject()
            ), 17
        );
}
</script>

The KML file could refresh automatically in Google Earth but not in OpenLayers. That is the real problem.

And one huge complexion is I cannot visualize the KML layer in OpenLayers. Either I use static or passive KML file is quite successful. Anybody?

Thanks in advance.


texts: Correcting inaccurate position coordinates

texts:
How to project the journey onto the street?

The black polyline is inaccurate GPS data received from a mobile device, which positions the user trajectory inside the buildings, when it was actually on the road.
The gray square pin is a sample inaccurate position at 59.436543,24.742133
The green circled pin is the origin at 59.4367469,24.742495
The green starred pin is the destination at 59.436377,24.742235
The red crossed pin is what I would like to calculate as the more accurate position estimate?

here's the geojson data pictured above
What's the formula?
One way to get X would be to draw a straight line between the origin and the destination (the green line) and then draw a perpendicular line to the inaccurate position from that green line.. but I have no idea what algorithm is needed to implement that?
function align(inaccurLat, inaccurLon, originLat, originLon, destinLat, destinLon) {
  // ???
  return [alignedLat, alignedLon];
}

What code would return the coordinates of X using the known coordinates of the circled pin, the starred pin and the square pin?
The calculation should work regardless of where the inaccurate square position is: it might be on the left of the green line, but it might also be on the right of the green line, or in fact above the green line (above the origin).
The function would ideally be in JavaScript, but I could do the translation to JavaScript myself if needed.
It should work offline.

Looking for a simple map matching algorithm
Furthermore, I can draw many green lines, because I have many origin and destination pairs available to me as the reference points representing road segments that do not curve too much - as depicted in the image below:


texts: How to draw a point at an exact location with spplot?

texts:
I imported a Shapefile holding a quite detailed world map from: http://www.diva-gis.org

library(sp)
library(rgdal)
world = readOGR('./countries_shp/', 'countries')

Using R with the sp-Packages, plotting is really easy:

spplot(world, 'UNREG2')
plot(world)

Now I want to

draw a customizable grid
add points derived from coordinates

For example importing coordinates of Vienna (48°12′32″N 16°22′21″E) can be done using char2dms of the sp-package:

vienna.N = char2dms("48d12'32\"N")
vienna.E = char2dms("16d22'21\"E")

But how can a draw a point at this exact location on my map?

Also adding a grid is straight forward:

plot(gridlines(world), add=TRUE)

assuming, that the last plot was created with "plot(world)".
How can this be achieved with the spplot output?


texts: What is the easternmost point in the US?

texts:
In his blog for The Atlantic magazine, James Fallows claimed that Eastport, ME was the easternmost point in the US. A writer commented that no, the easternmost point in the US was somewhere in Alaska, as Alaska crosses the 180th meridian. Another correspondent (me) wrote that this was silly - that the easternmost point of a landmass is the extreme point of the landmass that can be reached by traveling in an easterly direction, or more precisely, the last point on the country found by sweeping a line of longitude in an easterly direction. (This is more complicated for countries that don't have contiguous territory, but for now let's assume the territory is contiguous (and yes, I know that Alaska borders Canada and not another US state)).

Is there any convention in the geography community that establishes the meaning of "easternmost", preferably one that can be cited?


texts: Change main properties of map composer from cache to render from Python Console

texts:
I'm trying to print out a map from Python console. The issue is that when I add a layer to the map composer, the property of the map is set to cache. However, I want to change it to render so that the map would look less blurry.

Below is my code to add a layer to the map composer:

map_settings = iface.mapCanvas().mapSettings()
c = QgsComposition(map_settings)

view = iface.createNewComposer()
view.setComposition(c)

composerMap = QgsComposerMap(c, 10 ,10, 160, 20)
composerMap.setBackgroundEnabled(True)

rect = QgsRectangle(map_settings.fullExtent())
composerMap.setNewExtent(rect)
composerMap.updateItem()
c.addComposerMap(composerMap)

I am wondering if there is any way to set the main property of map as render instead of cache from Python console.


texts: Is there a way to export OpenStreetMap public transport routes to a KML or KMZ file? I'm trying to view the routes in another application (Google Earth) which supports KML data, but doesn't support the OpenStreetMap API.


texts: Export as Python script from ModelBuilder never completes?

texts:
The script below was exported from ModelBuilder and runs within ModelBuilder in about 7.5 seconds.

However when I export it and try to run it as a standalone script in either PyScripter or the default Python IDLE it creates a schema lock on the DB and never ends...I have let it run for an hour.

There is no error it simply never completes the task.

Any ideas on how to trouble shoot this or why it works in ModelBuilder but not as a standalone?

I'm using ArcGIS 10.2 for Desktop with a Basic level license.

    import arcpy

    arcpy.env.overwriteOutput = True

    # Local variables:
    Wells_csv = "c:\\xx\\Wells.csv"
    Wells_Layer = "Wells_Layer"
    wells = "c:\\xx\\wells"

    # Process: Make XY Event Layer
    arcpy.MakeXYEventLayer_management(Wells_csv, "Surface_Longitude", "Surface_Longitude", Wells_Layer, "GEOGCS['GCS_North_American_1927',DATUM['D_North_American_1927',SPHEROID['Clarke_1866',6378206.4,294.9786982]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98305509728916E-09;0.001;0.001;IsHighPrecision", "")

    # Process: Copy Features
    arcpy.CopyFeatures_management(Wells_Layer, wells, "", "0", "0", "0")


texts: LiDAR data download, merge and mosaic (Minnesota state)

texts:
I have been using the toolbox found here to download LiDAR data for MN (Minnesota-US).

The tool is very convenient for downloading multiple tiles but the merge/mosaic function of the tool is failing. The result is all the tiles are stored individually in separate gdb's which makes it very time consuming to extract in order to mosaic them on my own.

Is the merge/mosaic function of the tool failing because the tool is designed for 10.1 and I'm using 10.2? Anyone else have this problem or a solution?
Is there a way to combine/extract all of these tiles from multiple gdb's into one gdb in order to mosaic all the tiles?


texts: How to know if a point is on linestring or near one (<0.01units) using GeoPandas

texts:
I want to know if a point intersects/touches/on the line/linestring or if the point is less than 0.01 mts (units) from linestring.
import geopandas as gpd
import matplotlib.pyplot as plt

lines_gdf = gpd.read_file('line.shp')
points_gdf =  gpd.read_file('points.shp')

pip_mask = points_gdf.intersects(lines_gdf.loc[0, 'geometry'])

post running the script I can see all the data in pip_mask is false.
So if intersects/touches doest work is there any way to check if the point is <0.01 mts form the line?


texts: Load WFS from ArcGIS Online on OpenLayers 3

texts:
My goal is to load several shapefiles and display them using OpenLayers 3 as the basemap. I now have these shapefiles as a Web Feature Service (WFS) on ArcGIS Online. Below is what I have done (starting with one WFS).

// arcgis start
  var geojsonFormat = new ol.format.GeoJSON();

  var vectorSource = new ol.source.Vector({
    loader: function(extent, resolution, projection) {
      var url = '// https://dservices1.arcgis.com/DbPykcCwUUYq5zKg/arcgis/services/urban/WFSServer?service=wfs&request=getcapabilities';
      $.ajax({url: url, dataType: 'jsonp', jsonp: true});
    },
    strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
      maxZoom: 19
    }))
  });
  window.loadFeatures = function(response) {
    vectorSource.addFeatures(geojsonFormat.readFeatures(response));
  };

  var centerlocations = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style()
  });

  // arcgis end

My browser console shows ol.tilegrid.createXYZ is not a function. How do I go about successfully loading my shapefiles?


texts: What are definitions of watershed segmentation, local maxima, variable window search, regarding LiDAR for vegetation?

texts:
I am not familiar about Lidar. I want to measure trees height and number of trees.

I guess, watershed = to locate the trees, local maxima = to measure trees height. I used these methods in SAGA GIS.

I want to compare with Lidar Analyst. What are methods in Lidar analyst so that it can automatically count trees height and trees position?


texts: Construction of a square grid from a given centre (lat,lon) and direction of the sides

texts:
Let us consider a point (longitude, latitude) on Earth (location).

How to construct a square grid around this point (taking this point as the centre of the square or intercept of the two diagonals) with size of the square as 0.1 degree times 0.1 degree and direction of the one side of square is theta degree?

For example, let the sample point to be (88.38045, 21.04247) and theta = 70 degree. Here is a rough plot for clarification.

Note: This direction is computed from the 84 E meridian.


texts: Slow performance of ST_CONTAINS on MySQL 5.7

texts:
I'm using MySQL 5.7 and I have two tables, one containing the geographies of every census block group in the US, for a total of 223,325 rows.

I have another table containing approximately 400 custom geographies. For a given custom geography, I want to retrieve a list of all census block groups that reside within it.

However, the query seems to be taking far too long (well over 30 minutes before I had to cancel it)

I'm using the following query:

SELECT
    geo_bg.OGR_FID
FROM
    geo_custom
        JOIN
    geo_bg ON st_contains(geo_custom.shape, geo_bg.shape)
WHERE
    geo_custom.OGR_FID = 213;

And the output from EXPLAIN is

+----+-------------+------------+------------+-------+---------------+---------+---------+-------+------+----------+-------------+
| id | select_type | table      | partitions | type  | possible_keys | key     | key_len | ref   | rows | filtered | Extra       |
+----+-------------+------------+------------+-------+---------------+---------+---------+-------+------+----------+-------------+
|  1 | SIMPLE      | geo_custom | NULL       | const | OGR_FID,SHAPE | OGR_FID | 4       | const |    1 |   100.00 | NULL        |
|  1 | SIMPLE      | geo_bg     | NULL       | range | SHAPE         | SHAPE   | 34      | NULL  | 9952 |   100.00 | Using where |
+----+-------------+------------+------------+-------+---------------+---------+---------+-------+------+----------+-------------+

I have a spatial key on the shape column in both geo_bg and geo_custom

Is there something obvious I'm missing, or is this type of query very expensive and going to take a long time no matter what?


texts: Creating file geodatabase and using it as variable to create file path?

texts:
It's been a while since I've done arcpy scripting so it feels kind of new again. I'm working a script that is asking for these 5 parameters:

workspace = arcpy.GetParameterAsText(0)
RouteLayer = arcpy.GetParameterAsText(1)
StateLayer = arcpy.GetParameterAsText(2)
StateTaxDistricts = arcpy.GetParameterAsText(3)
NewGDB = arcpy.GetParameterAsText(4)

arcpy.CreateFileGDB_management(workspace, NewGDB)

What I want to do is to use the new GDB as a variable and create a location for outputs. I want to do something like this:

newLocation = arcpy.CreateFileGDB_management(workspace, NewGeodatabase)
clipLayer = arcpy.Clip_analysis(RouteLayer, StateLayer, newLocation + name + "Clip")

I've got a few other tools that are using the same format for the output. I have a clip, intersect, split line at point, and spatial join.

When I run my tool, it crashes and calls out this line: newLocation = arcpy.CreateFileGDB_management(workspace, NewGeodatabase)

What am I doing wrong?


texts: Maintaining reference to feature layer (rather than its name string) throughout Python script?

texts:
When creating a Layer in arcpy, the syntax uses a string to denote the layer name, as shown in the Make Feature Layer sample script:

import arcpy

arcpy.env.workspace = "C:\Program Files (x86)\ArcGIS\Desktop10.2\TemplateData\TemplateData.gdb"
arcpy.MakeFeatureLayer_management("city", "citiesLyr")

print(type("citiesLyr")) # => string

As far as Python is concerned, "citiesLyr" is just a string (right?). Can I instead define the feature layer as a variable?

This is an attempt to simplify a more complicated scenario. I'm trying to use the feature layer in a script which involves multiprocessing, loops and functions, and I'm finding that the reference to the layer is being lost, and all I'm left with is a string called "citiesLyr". How can I keep the reference to the feature layer throughout my script?

The example below doesn't do anything - the point is to test how to gain access to the feature layer from within the function, which has been called within the multiprocessing environment:

import multiprocessing, arcpy

def doCity(lyr):
    #How to get access to the cities layer from this function?
    print(lyr)
    desc = arcpy.Describe(lyr) # <= this fails because the layer is just a string
    print(desc.SpatialReference)

if __name__ == '__main__':

    #Create a feature layer from the cities
    arcpy.MakeFeatureLayer_management("C:\Program Files (x86)\ArcGIS\Desktop10.2\TemplateData\TemplateData.gdb\city", "citiesLyr")

    #Create a list, so we can use the multiprocessing function
    citiesList = []
    for i in range(0,5):
        citiesList.append("citiesLyr")

    # Create a pool class and run the jobs
    pool = multiprocessing.Pool()
    pool.map(doCity, citiesList)

    # Synchronize the main process with the job processes to ensure proper cleanup.
    pool.close()
    pool.join()

(The cities list itself is also pointless - it's only there to use  multiprocessing)


texts: Select the smallest value from a selection using ArcMap SQL selection

texts:
I am trying to select the smallest Directionality value based on a selection of two Workingspace values:

In the above picture, there are two "Jeffery's" selected, I want to then select the Jeffery with Directionality 2 (I.E. the smaller one) without using = 2. I know there is a max/ min SQL expression:

DIRECTIONALITY = (SELECT MIN( DIRECTIONALITY) FROM TableName)

But this returns no values when choosing the select from current selection method.
HOWEVER, if I create a new selection, this SQL does return a selection of 1 (it selects the row where Directionality is 1).

Is there a way to tell ArcMap to select the smallest value of a previous selection based on a different field? Specifically in Model Builder.


texts:
I need to create a chart from the result of a query task using Esri JS API.  Note, the chart does not have to be created with Esri API, but the query runs in that.  The result is written to the DOM and can be displayed using an .innerHTML.

So I basically have a number (Query Result) that is given to the DOM that I would like to use to build a chart.  The fun thing is that its a time aware DynamicService and this makes great use of the
timeslider widget to create a great visualization of development over time.  I now wish to add to this with a statistical review.  I am able to return a total built number to the DOM and display that in a DIV, as shown in the JSFIDDLE.  So its a cool visualization and provides some stats.  the query is running on the same same source data, but on a REST that is not M-enabled.

What I am looking for is to provide a line chart that adds records to it as the number increases - query results returned to the DOM.  I have no idea where to start.


texts: Creating style rules based on more than one criteria

texts:
I have one shapefile containing parts of hidrographic basins on the brazilian state of Espírito Santo. I want to generate a classification style/label based on two attributes:

1st: The percentage of forest on the landscape of each basin. The data should be divided between 0 - 10% and 10-50% landscape coverage (there's no more than 50% coverage on this shapefile)

2nd: The Patch Connectivity index generated from my analysis, ranked in crescent order.

Basically, the shape should be divided based on the percentage of forestation and each class would recieve a different color gradient, that would rank the PC index. From what I know about QGIS, this should be done using a rule based styling. I tried it using the percentage of forest as my first classifier and then adding categories to it. Here's what I came up with (let's try to overcome the language barrier here):

Problem is, after I apply these settings, the shapefile appear as if it's only being classified using the percentage of landcape.

I have tried deactivating the plain color label, but the whole shapefile disappears. If I change the style of the 1st level label to a translucid one, there's nothing beneath. I also tried changing the symbol levels (the order they're drawn) and nothing happened. Any thoughts on how to solve this?


texts: How to recreate past landscapes on ArcGIS 10.1?

texts: I am trying to recreate the landscape as it may have been around 4000 BC. I have an idea about what tree species existed in particular regions of Scotland (by reading pollen studies). Now I need to somehow add trees to the Land-Form PANORAMA raster (taken from Ordnance survey UK) to visualize the landscape. I am aware of adding layers (with different datasets - for instance geological information) but I am not sure how I can add tree cover in ArcGIS as I don't have this information in the form of a layer which can be added to the raster.I was wondering whether anyone could point me in the right direction. Any help or advice would be greatly appreciated.


texts:
I'm using QGIS 3.0 with the expanded SAGA options. I created a training file in SAGA standalone that I am using, but would prefer the GRASS kappa comparasion options. I have already performed the classifications in SAGA standalone but want to use qgis.
I'm just a little confused about if I'm going to perform the classification and select every band individually and if that will give me 6 outputs. It's a landsat 8 image, 30 km by 50km.
If anyone knows of additional documention, tutorials about the saga toolbox within QGIS.


texts: Polygons from two different shapefiles that share common border in QGIS?

texts:
I have two shapefiles, which are village administrative units in this case.

Given that these two shapefiles share a common boundary/border is it possible to get the border units (villages) on both side of the border?


texts: Prevent Export/Save As of Vector Layers/Shapefiles

texts:
I have created a database using postgreSQL and created multiple users assigned with different roles

Some of the users can only select and view the data and not update or delete any record. It works perfectly both in postgres and QGIS environment. However, the user can easily export the table into a shapefile in QGIS and then edit or update whatever he wants to do with the exported shapefile.

Wondering if there is a way I can lock the shapefile/table so that user cant save or export the postgis table to a shapefile in QGIS? I read something about disabling the SAVEAS button in QGIS Registry, but is this the only possible way?


texts: QGIS turn right click -- "Toggle Editing" off

texts:
In QGIS3, Interface customization I have turned off the mEditMenu but when a user right clicks on a dateset they still see "Toggle Editing".

Is it possible to remove this? and also any keyboard shortcuts that allow the same.


texts: Python Floating Point Error when using SCP plug-in (QGIS 3.6)

texts:
I receive a python error every time i try to save a temporary ROI to training input in the Semi-automatic classification plug-in in QGIS 3.6.  I'm unable to add any training ROI to training input.

python error below

i've tried:

uninstalling, and reinstalling the plugin
installing a new version of qgis (was originally in 3.4)
used other images

Any one have a similar problem and found a solution?

Error Message below:

2019-05-21T19:02:02     WARNING    Traceback (most recent call last):
          File "C:/Users/ncu116/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\dock\scpdock.py", line 2503, in saveROItoShapefile
          cfg.utls.calculateSignature(cfg.shpLay, cfg.bandSetsList[cfg.bndSetNumber][8], [self.ROILastID], cfg.ROIMacroID, cfg.ROIMacroClassInfo, cfg.ROIID, cfg.ROIInfo, 50, 40, "No", "No", UID)
          File "C:/Users/ncu116/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\core\utils.py", line 1739, in calculateSignature
          ROIsize = cfg.utls.getROISize(bX, 1, rStat[0], rStat[1])
          File "C:/Users/ncu116/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\SemiAutomaticClassificationPlugin\core\utils.py", line 5672, in getROISize
          s = ( (a >= min) & (a              FloatingPointError: invalid value encountered in greater_equal


texts:
I'm trying to use a publicly available data set (the Australian National Vegetation Information System) in QGIS, but it has been set up for use in ArcGIS (.lyr file).

I tried to convert the file using http://wald.intevation.org/projects/arcmap2sld/ but kept getting an error saying it was an unsupported file type.

I can open the individual data layers directly in QGIS, however, there doesn't seem to be any coherency to the way the data layers are named, and hence I'm not sure what the differences between them are, or whether I'm using the appropriate one for my project.

Can anyone recommend another approach for using the data in QGIS and/or advise on how to interpret the data package?

Before asking I reviewed:

this article on the difference between .lyr files and .qlr files,
and
this question about converting .lyr files

I also looked at the slyr tool, but couldn't figure out how to install it (sorry, newbie here!)


texts:
I am working with global applications and decided to use the Python Earth Engine API to collect some data I will need. I am exporting annual, global NDVI composites based on Landsat. Since this involves a lot of data, I split the task in small tiles leading to about 1000 downloads to be conducted. Due to the high numbers of tasks, I run into a 500 error after a while, which is becoming quite annoying.

My question is this: after downloading a file, is there a way to clear the memory used to derive that file?


texts: Setting field visibility in ArcMap using C#?

texts:
I want to set certain fields of a featurelayer to be visible on startup. how would I go about changing the field visibility of certain fields?

if (featureLayer.Name == "LayerName")
{
    // Code to set certain fields to visible
}

I tried using the IFields

IFields fields = featClass.Fields;
IField field = null;
for (int i = 0; i < fields.FieldCount; i++)
{
    field = fields.get_Field(i);

    if (field.Name == "ObjectID")
        // field.Visible = false;
}

but I dont see a visibility property

UPDATE: Here is the equivalent code based on Hornbydd's answer below in C#

void SetFieldVis()
{
    // Get Document
    IMxDocument pMXD = ArcMap.Application.Document as IMxDocument;

    // Get Map
    IMap pMap = ArcMap.Document.FocusMap;

    // Get first layer in TOC -> in my case set up an enum and loop through the layers
    UIDClass uid = new UIDClass();
    uid.Value = "{40A9E885-5533-11d0-98BE-00805F7CED21}";
    IEnumLayer enumLayers = map.get_Layers(uid, true);
    IFeatureLayer pFeatureLayer = enumLayers.Next() as IFeatureLayer;

    while (pFeatureLayer != null)
    {
        if (pFeatureLayer.Name == "LayerName")
        {
            // Get a pointer to the fields
            ILayerFields pLayerFields = featureLayer as ILayerFields;

            // Get the fourth field by it's index number and set its visibility to False
            pLayerFields.FieldInfo[3].Visible = false;
        }
    }
}


texts: OSM - Extract Streets from OSM Data (With/Without Python)

texts:
I'm currently in the works of creating a pythonic way of extracting street names from an OSM dump for the Province of Quebec. Originally, I wanted to do this in Python to make it easier to integrate the process in future code, but I'm failing miserably, and after weeks of on and off research without getting anywhere, so now I'll honestly take what I can get.

What's the best way to extract all the street names from a dump such as these?
http://download.geofabrik.de/north-america/canada/quebec.html

Here's my original code if anyone can make it work. At the moment this code maxs my 16GB of ram and 32GB of swap in about a minute before it evens begins going through each element.

from lxml import etree
import xmltodict, sys, gc
from pymongo import MongoClient

#Ultimate fix of life.
reload(sys)
sys.setdefaultencoding("utf-8")

client = MongoClient()
db = client.re
streetsDB = db.streets

hwTypes = ['motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'pedestrian', 'unclassified', 'service']

#Enable Garbadge Collection
gc.enable()

def fast_iter(context, func):
    print 'Begun processing'
    placement = 1
    # http://www.ibm.com/developerworks/xml/library/x-hiperfparse/
    # Author: Liza Daly
    for event, elem in context:
        placement += 1
        print placement
        func(elem)
        elem.clear()
        while elem.getprevious() is not None:
            del elem.getparent()[0]
    del context

def process_element(elem):
    data = etree.tostring(elem)
    data = xmltodict.parse(data)
    keys = data['way'].keys()
    if 'tag' in keys:
        if isinstance(data['way']['tag'], dict):
                    if data['way']['tag']['@k'] == 'highway':
                        if data['way']['tag']['@v'] in hwTypes:
                            streetsDB.insert(data)
        else:
            for y in data['way']['tag']:
                if y['@k'] == 'highway':
                    if y['@v'] in hwTypes:
                        streetsDB.insert(data)
                        break
    del data
    del keys
    gc.collect()

context = etree.iterparse('quebec-latest.osm', tag='way' )
fast_iter(context,process_element)


texts:
We have the QWC2 client installed on our server and it opens correctly in the client with the standard Swiss layers.
We have a QGIS WMS Server running with a number of layers.
We want to know how to configure the themes.json so that it serves our layers.
The documentation says:

themes.json - Full theme configuration, autogenerated from themesConfig.json

Creating themesConfig.json looks laborious given we have a large number of layers.
Is there any sort of utility that can generate a basic themesConfig.json from a GetCapabilities XML generated by the WMS server?


texts: Find and remove portions of polygons from one geometry column that overlap from another geometry column in PostGIS

texts:
Let's say I have 2 VIEWS. Each has a geometry column with polygons that overlap, approximately like so:

I'd like to delete the portions of of the green polygons where there is any overlap with the grey ones, approximately like so:

Both VIEWS have polygons that extend out of the boundaries of the example images to the left and right, respectively.

I took a shot at this with the following query (based on a similar example from a book):

SELECT my_value, ST_Intersects(g1.geom1, g1.geom2) As they_intersect
, ST_Difference(g1.geom1, g1.geom2) As difference_polygon
FROM (SELECT
view_green.geom As geom1,
view_grey.geom As geom2,
view_green.my_value As my_value
FROM view_green, view_grey) AS g1;

(keeping a value associated with each polygon as well)

The output from that is the same as the VIEW and does not work (I'm using openlayers to visualize at the VIEWS/results). So as a test, I wrote a query to return only the overlapping polygons and delete the areas that overlap:

SELECT
view_green.my_value,
view_green.timestamp,
CASE WHEN (view_green.timestamp < view_grey.timestamp)
     THEN ST_Difference(view_green.geom, view_grey.geom)
END
FROM
view_green,
view_grey
WHERE ST_Intersects(view_grey.geom, view_green.geom);

(You'll notice the .timestamp because I'm trying to delete overlap from a view that is 'older'. Eventually I'd like to have a TRIGGER similar to this solution and skip the timestamp comparison)

The result of that query does not work fully but performs the deletion correctly (with st_difference()) on only the first polygon, approximately like so:

Where am I going wrong?


texts: Join field in modelbuilder won't work, but will manually

texts:
I'm trying to do a 'join field' in a model. When I do the join simply not in a model, no problem.  When I try to use the 'join field tool' in a model, I get the ""<sheet> does not have OIDs.

I have a very simple feature class in a GDB, trying to join with an Excel sheet. No problem to join. Only in the model.

Anyone got a clue?

Both files comply with the 'join rules' in file names etc.


texts:
I have registered the postgresql database with arcgis desktop 10.1.

In ArcMap, When I try to drag and drop one of the tables to the map I get this error:

Could not add the specified data object to the map

I want to Add data as a layer to the map using the data tables in postgresql database rather than using add layer plugin and import .csv file form another location.

How can I solve this?


texts:
Hi I like to build a GUI for map like this one: https://www.youtube.com/watch?v=uzuiEnbVU2U
but using OSM data.
I know the syntax of Tk, but have a problem to put the map inside. Need some advices/tips to load the map tile inside the window asynchronous and make it scrollable.


texts:
I'm using d3.js to display isobars from GeoTIFF in Leaflet.
The problem I would like to solve is the smoothing of those isolines.
Here is my code:
var xhr = new XMLHttpRequest();
    xhr.open('GET', `./20210129T030000Z_1002.000.tif`, true);
    console.log(xhr)
    xhr.responseType = 'arraybuffer';
    xhr.onload = function () {
        console.log(this.response)
        var tiff = GeoTIFF.parse(this.response);
        var image = tiff.getImage();
        var tiffWidth = image.getWidth();
        var tiffHeight = image.getHeight();
        var rasters = image.readRasters();
        var tiepoint = image.getTiePoints()[0];
        var pixelScale = image.getFileDirectory().ModelPixelScale;
        var geoTransform = [tiepoint.x, pixelScale[0], 0, tiepoint.y, 0, -1 * pixelScale[1]];
        console.log(geoTransform)

        var pressData = new Array(tiffHeight);
        for (var j = 0; j < tiffHeight; j++) {
            pressData[j] = new Array(tiffWidth);
            for (var i = 0; i < tiffWidth; i++) {
                pressData[j][i] = (rasters[0][i + j * tiffWidth] / 100).toFixed(2);
            }
        }
        var intervalsPress = [970, 972, 974, 976, 978, 980, 982, 984, 986, 988, 990, 992, 994, 996, 998,
            1000, 1002, 1004, 1006, 1008, 1010, 1012, 1014, 1016, 1018, 1020, 1022, 1024, 1026, 1028];
        var isobars = rastertools.isolines(pressData, geoTransform, intervalsPress);
        console.log(isobars)

        isobarsLayer = L.geoJSON(isobars, {
            style: {
                "color": "#f00",
                "weight": 1,
                "opacity": 0.65
            }
        });
        isobarsLayer.addTo(map)
    }
    xhr.send()

Here is the output:

Can I use any interpolation like bilinear or bicubic? If yes, how?
Or more easy, can I smooth the lines?


texts: QGIS does not show FID field of vector data

texts: Shapefiles (vector layers) has FID or OID fields associated with each row (feature). This information is readily available in the attribute table viewed in Arcmap. However, QGIS Desktop does not seem show the FID field. Does anyone know how to view FID field in QGIS Desktop?


texts: Calculating the mean NDVI value between top 20% to 30% lowest NDVI in the geomerty for all imagecollection

texts:
The purpose is to calculate the representative NDVI value in the geometry for all imagecollection.
The representative NDVI in one image should satisfy the requirements that

To ascend the NDVI in one image

To extract the NDVI between top 20% and top 30% from STEP ONE

To calculate the mean NDVI from STEP TWO.

To add the mean to the imagecollection.
    geometry =
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-117.77231150781333, 40.25497601925202],
          [-117.77231150781333, 39.84710661565274],
          [-117.08017283593833, 39.84710661565274],
          [-117.08017283593833, 40.25497601925202]]], null, false);

var imageCollection = imageCollection   .filterDate('2014-01-01', '2015-01-01')   .filterBounds(geometry)

Map.addLayer(imageCollection,null,'original imagecollection')

var addNDVI = function(image) {   var ndvi = image.normalizedDifference(['B5','B4']).rename('NDVI').float(); // first the NIR band, then the red band   return image.addBands(ndvi); };

var imageCollection = imageCollection.map(addNDVI)

print('added NDVI imagecollection',imageCollection )```

code


texts: Importing CSV as a point layer via the Python console. No point layer is loaded. How can I fix my code?

texts:
Code
import os
from qgis.core import (
    QgsVectorLayer
)

uri = r"file:\\\C:\Users\imper\Desktop\XYbuffer_tool\Calapan_XYpoints_comma?delimiter={}&crs=epsg:32651&xField={}&yField={}".format(",", "x_long", "y_lat")
vlayer = QgsVectorLayer(uri, "test", "delimited text")

QgsProject.instance().addMapLayer(vlayer)

From the PyQGIS documentaiton https://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/loadlayer.html#vector-layers
Result
No point layer is loaded
exec(open('C:/Users/imper/AppData/Local/Temp/tmpo4j1vq9h.py'.encode('utf-8')).read())

CSV file

name,corner,x_long,y_lat Shell Calapan,1,121°12'07.5"E,13°25'42.6"N
Shell Calapan,2,121°12'09.1"E,13°25'38.8"N  Shell
Calapan,3,121°12'06.0"E,13°25'37.6"N Shell
Calapan,4,121°12'04.3"E,13°25'41.5"N


texts: Postgis 2.0 upgrade problem with postgres data?

texts:
I have a serious problem:
on Wednesday, April 2 OpenSUSE 12.1 is updated, and the package was to update the PostGIS-a version 1.5.3 to 2.0.0.
  After this operation, Postgres does not display records.
  In pgAdmin count the records, but seeing them in a table view, and by SELECT is a problem.
  I can not do dump:
pg_dump-h 127.0.0.1-p 5432-U postgres-Fc-b-v-f "/ home / admin / olddb.backup" Tomice

I get the message:
  pg_dump: SQL command failed
  pg_dump: Communication error with server :: ERROR: can not access the file "$ libdir/postgis-1.5 ': No such file or directory
  pg_dump: Your command: COPY public.armatura (gid, type, sub_typ, description, wys_g, wys_d, operator, the_geom) TO stdout;
  pg_dump: *** stopped because of an error

In view of the tables in PgAdminie my tables, I have a good number of records, but when you call the data is still the same.

I noticed that almost all functions have the following address:
  "$ libdir/postgis-1.5" that during the upgrade-u has been replaced by / PostGIS-2.0.
  I could perhaps remove the column geometry, but at the base of my ...

I attach screenshots of pgAdmin.
  Command did not help, sql for good-sounded names: postgis_upgrade_2.0.sql etc. ..
  Suddenly all the repositories of RPMs flown for PostGIS 1.5.3, and I can not go back to the old version ...
  At stake is 1 working day to 15 people. If anyone has any idea - I'll be bound.

I have no idea to solve this problem.


texts: Calculating distances between selected points from a single point in another layer in QGIS

texts:
I'd like to calculate the distances for all competitors to a proposed operation within a certain radius.
Layers are:

Competitors
Possible Competitors
Subject Property
Market Radius

Ultimately, I want to be able to add a "weight" to each competitor based on distance from the subject property (red spot in the center of the radius).

Any thoughts?

I'm working in QGIS.


texts:
I have recently been looking at airlines web pages that display their routes departing a certain city to all other cities they service.  I would like to be able to create similar curved routes between points.  Has anyone created scripts or functions that will generate the curved arcs like those displayed in this example?

In PostGIS, is there an implementation of the ST_MakeLine that would allow you to specify the amount of curve to use when connecting 2 points?

While I am currently using PostGIS and QGIS, I would welcome hearing about other software options that might be able to create the same appearance.


texts:
I am trying to create a map which shows points that are clustered (or not depending on the distance between the points).

So, I have the following query:

  SELECT
  ST_NumGeometries(gc) as num_geometries,
  ST_AsText(gc) AS geom_collection
  FROM (
    SELECT unnest(ST_ClusterWithin(coordinates, 0.0035)) AS gc FROM (
        SELECT DISTINCT ON (coordinates) coordinates FROM chargepoints
            WHERE coordinates && ST_MakeEnvelope(4.281578063964845, 51.81646817257424, 4.635200500488282, 51.99460541491795, 4326)) AS points
) as result;

I am using ST_MakeEnvelope to retrieve charge points within that rectangular map.
ST_ClusterWithin makes geometry collections where each collection holds a set of geometries:

2   GEOMETRYCOLLECTION(POINT(4.29696 51.85363),POINT(4.29721 51.85387))
12  GEOMETRYCOLLECTION(POINT(4.303751 51.854067),POINT(4.30149548 51.8520582),POINT(4.302743 51.853874),POINT(4.303193 51.854307),POINT(4.303058 51.854648),POINT(4.30405397 51.8524191),POINT(4.3007711 51.8529992),POINT(4.30717715 51.85126392),POINT(4.30502108 51.85145128),POINT(4.30586433410644 51.8506546020508),POINT(4.305203 51.852615),POINT(4.305643 51.851641))
1   GEOMETRYCOLLECTION(POINT(4.306028 51.836928))

What can or should I do with each collection to eventually get circled clusters on an map?

Please let me know if you need clarifications.

Using type geometry and SRID 4326, coordinates are in degrees (lat lon)


texts: Is it possible for Orfeo toolbox to do auto object identification on a 3band RGB raster image?

texts:
My intention of using Orfeo Toolbox is to find all the open air basketball courts on a raster image with 3 bands RGB.

Does OTB require a multi spectral band image or is 3 band image good enough to perform auto object identification.


texts: (QGIS) Shapefile with attribute table: Convert it to raster burning values of each polygon to the raster. "Heatmap"

texts:
I am working on QGis with a shapefile with an attribute table, I have created this shapefile through "Vector->Analysis tools->Count Points in Polygon", so this new shapefile got the count in the attribute table.

I want to convert now this shapefile to a "heatmap" in raster format.

I am trying to do it with "Raster->Conversion->Rasterize", burning the values of the attribute table

In those polygons with number=0 QGIS doesn't draw the polygons in the raster.

I would like to get all those polygons filled in the raster, I know I could erase all those polygons in the attribute table with data, then change all the 0s to 1, and then rasterize it, merging both rasters (the original with empty polygons and the new one with only the "empty" polygons filled).
I think this is just too tricky .

Is there any other solution?


texts: What could replace the Traversability Index (TRAVERSE) from previous Fragstats versions?

texts:
I need to analyse the connectivity of wetlands, I believe that the Traversability Index from Fragstats is the one that best describes my intentions because it considers not only the ecological cost of different landcovers but also the limit of an ecological budget. However this metric no longer exists in the program. Is there another metric within Fragstats that could perform the same task as in the 4.2 version of the program or is there any other GIS tool?

In order to explain my situation, here is how the metric is described

  This metric uses a resistance-weighted spread algorithm to determine the area that can be reached from each cell in a focal patch. The focal cell gets an organism-specific “bank account,” which represents, say, an energy budget available to the organism for dispersal from the focal cell. The size of the account is selected to reflect the organism’s dispersal or movement ability, and translates into ecological neighborhood size around a focal cell. Each patch type (including the focal patch type) is assigned a cost, based on a user-specified resistance matrix. Specifically, relative to a each focal patch type, each patch type is assigned a resistance coefficient in the form of weights ranging from 1 (minimum resistance, usually associated with focal patch type) to any higher number that reflects the relative increase in resistance associated with each patch type. The index is computed at the individual cell level as follows. The metric is computed by simulating movement away from the focal cell in all directions, where there is a cost to move through every cell. Even a cell of the same patch type will have a cost, usually set to 1, so that under the best circumstances (i.e., minimum resistance), there will exist a maximum dispersal area based on the specified account or energy budget. Note that assigning a (small) cost for traveling through the focal community (typically a cost of 1) results in a linearly decaying function. Moving through more resistant cells costs more and drains the account faster. Thus, depending on the resistance of the actual landscape in the vicinity of the focal cell, there will be a certain area that a dispersing organism can access. This area represents the least-cost hull around the focal cell, or the maximum distance that can be moved from the cell in all directions until the “bank account” is depleted. This dispersal area, given as a percentage of the maximum dispersal area under conditions of minimum resistance, provides a measure of the traversability of the landscape in the vicinity of the focal cell. Averaging this index across cells at the patch, class, or landscape level provides an index of traversability

I would like to create something like this picture from FRANK, S. Development and Validation of a Landscape Metrics Based Approach for Standardized Landscape Assessment Considering Spatial Patterns. 2014.


texts: Accessing OGR driver options using OGR/GDAL in Python?

texts:
Using OGR in Python, I am attempting to modify the driver settings for S57. I am able to read the meta data from the driver as follows:

import ogr

driver = ogr.GetDriverByName("S57")
driver.GetMetadata()

Which returns all the meta-data for the driver, I have removed some of it for the sake of the post length.

{'DCAP_CREATE': 'YES',
 'DCAP_OPEN': 'YES',
 'DCAP_VECTOR': 'YES',
 'DCAP_VIRTUALIO': 'YES',
 'DMD_CREATIONOPTIONLIST': "<CreationOptionList>   <Option name='S57_EXPP' "
                           "type='int' description='Exchange purpose' "
                           "default='1'/>   <Option name='S57_INTU' type='int' "
                           "description='Intended usage' default='4'/>   "
                           "<Option name='S57_EDTN' type='string' "
                           "description='Edition number' default='2'/>   "
                           "<Option name='S57_UPDN' type='string' "
                           "description='Update number' default='0'/>   "
                           "<Option name='S57_UADT' type='string' "
                           "description='Update application date' "
                           "default='20030801'/>   <Option name='S57_ISDT' "
                           "type='string' description='Issue date' "
                           "default='20030801'/>   <Option name='S57_STED' "
                           "type='string' description='Edition number of S-57' "
                           "default='03.1'/>   <Option name='S57_AGEN' "
                           "type='int' description='Producing agency' "
                           "default='540'/>   <Option name='S57_COMT' "
                           "type='string' description='Comment' default=''/>   "
.........
'DMD_EXTENSION': '000',
 'DMD_EXTENSIONS': '000',
 'DMD_HELPTOPIC': 'drv_s57.html',
 'DMD_LONGNAME': 'IHO S-57 (ENC)',
 'DMD_OPENOPTIONLIST': "<OpenOptionList>  <Option name='UPDATES' "
                       "type='string-select' description='Should update files "
                       "be incorporated into the base data on the fly' "
                       "default='APPLY'>    <Value>APPLY</Value>    "
                       '<Value>IGNORE</Value>  </Option>  <Option '
                       "name='SPLIT_MULTIPOINT' type='boolean' "
                       "description='Should multipoint soundings be split into "
                       "many single point sounding features' default='NO'/>  "
                       "<Option name='ADD_SOUNDG_DEPTH' type='boolean' "
                       "description='Should a DEPTH attribute be added on "
                       "SOUNDG features and assign the depth of the sounding' "
                       "default='NO'/>  <Option name='RETURN_PRIMITIVES' "
                       "type='boolean' description='Should all the low level "
                       'geometry primitives be returned as special '
                       "IsolatedNode, ConnectedNode, Edge and Face layers' "
                       "default='NO'/>  <Option name='PRESERVE_EMPTY_NUMBERS' "
                       "type='boolean' description='If enabled, numeric "
                       'attributes assigned an empty string as a value will be '
                       "preserved as a special numeric value' default='NO'/>  "
                       "<Option name='LNAM_REFS' type='boolean' "
.......

Using the following functions

driver.GetMetadataItem("DCAP_OPEN")
driver.SetMetadata({"DCAP_OPEN" : "NO"})

I am able to obtain the "higher"-level options and set them accordingly. However, I'm interested in modifying options within the DMD_CREATIONOPTIONLIST and DMD_OPENOPTIONLIST, such as setting the UPDATES option in the latter. How do I check the status of said option and setting it afterwards?

Edit: Specifically the options on page 301 "S57 Control Options" link


texts: I'm working with MODIS AOD L2 data and would like to create daily averages and monthly averages from the data. I'm aware of the L3 product but it doesn't allow me to create an 18 day average of just weekend averages for example. Therefore, I'm wondering if anyone knows how to merge MODIS AOD L2 products to create L3 daily averages.


texts: Changing rotation for each data driven page dynamically?

texts:
I have a .mxd using data driven pages. In the .mxd there is a shape containing foto locations displayed with an arrow. At the top of the arrow there is a label containing the foto id visualized, the roation of the arrow is the direction of view of each photo.
The problem now is, that the data driven pages have different rotations, so the labels are not positioned right in each page. I am using the adjusting placement by rotation setting.
How can I dynamically change the rotation for each data driven page for each lable?


texts:
My SOE extension crashed. Now it does not appear in extension list, but does so in Capitabilities of a service. Services with it give error on run. When I try to add it anew to the server, I get
"An extension with name 'xxx' is already registered for service type 'MapServer'"
error.

ArcGIS Installation 10.1. To default location.

Note: Solutions I got on the internet:
      a. Restart ArcGIS Server from Task Manager( restarted many times everything that can be restarted -> no effect )
      b. Remove and re-install ArcMap Server ( can't do it now, may be as last resort)


texts: How can I relate the osm2po osm_source_id to data in OpenStreetMap?

texts:
I created and imported OpenStreetMap (OSM) data into Postgres using osm2po:

e.g. java -Xmx8g -jar osm2po-4.6.9/osm2po-core-4.6.9-signed.jar config=osm2po.config.eu_xc prefix=eu tileSize=x,c europe-latest.osm.bz2

and psql -d osm2po_eu_xc -q -f "eu_2po_4pgr.sql"

which creates the table eu_2po_4pgr and has datasets with columns "osm_source_id", "osm_target_id".

I created and imported the exact same osm data (using the exact same bz2 file) using osm2pgsql (0.69)

e.g. time osm2pgsql -smuc europe-latest.osm.bz2   -d osm_eu

which creates table like planet_osm_point which has a "column osm_id".

I can't find the "osm_source_id" in say planet_osm_point.

What I would like is to join the osm2po table with the original OSM Data, in order to have the complete information.

Is that possible or am I missing something fundamentally?


texts: QGIS Connection with MSSQL localdb not possible correctly

texts:
I installed MSSql Server Express with localDB. Everything looks ok. The connection with the LocalDB is nessesary and the tables and views are displayed in the QGIS MSSQL Dioalog. A primary key is defiened! See here:

... the Connection Dialog worked fine

... the QGIS MSSQL Table Dialog

The features are not displayed, and the attributetable is displayed as empty
The attribute table is always empty, if I connect a table with no geometry

Some notes and questions:

Using QGIS 2.10.1

MS SQL Server Express localdb with SP2

DSN Connection String for localDB with:"DRIVER=SQL Server Native Client 11.0;Integrated Security=True;"

Is something wrong with the Connection String?
It is not possible to connect with the Native Client?
Is somethin wrong with the tabledeffinition?
Why the column primary key in the QGIS Dialog is empty?

...or it is not necessary to connect to a localDb


texts: When creating an expression in QGIS, is there a way to reference the next record and previous record?

texts:
I have a set of points from a GPX file, and would like to create a heatmap from them, with the points weighted on the amount of time spent at that point.

I want to create a calculated field in QGIS, where duration will equal, either T-t where 'T' is the timestamp at the next point and 't' is the timestamp of the current point; or even better, (T-t)/2 + (t-s)/2 where 's' is the timestamp of the previous point.

However I'm unclear on how to bring values from other records into an expression.


texts: polygon/polyline digitizing techniques for geological data

texts:
My current work involves extensive digitising of historic geology data - essentially the creation of new polygons and polylines based on overlain raster and vector layers. This is all done manually using QGIS 2.14.3. I was hoping to get some tips from the community to make this process more efficient and to attain a better outcome.

As it stands now, I am simply tracing features and creating polygons and lines by hand. I've found this quite a frustrating process in QGIS. How do I create perfectly snapped adjacent polygons? Is the only way to play with snapping options (tolerances)? I've tried this but it still seems a bit clunky when you have a large density of vertices, and I still seem to miss some, and/or snap to polygons close by that I want to keep separate. In MapInfo, I could roughly digitise the boundary of the neighbouring polygon (making sure it overlaps) and then erase this overlap with a click of the mouse.

Similarly, I cannot find a tool that will allow me to easily edit the shape of a polygon. In MapInfo, I could simply draw another overlapping (overlaying) polygon, and use the erase tool to erase its extent from the underlying polygon. In QGIS, all I've found that comes close is the add ring feature, however it must be totally enclosed in the polygon being edited. As a result, when I want to edit the shape of a polygon, I find myself manually selecting, moving and adding nodes. This takes forever in some instances, and you have to add and move many many nodes in order to get a smooth result. I probably could draw an overlapping polygon in another layer and then use the layer based geoprocessing tools (Difference?) but I'm hoping there is a more simple solution.

Can anyone help me? I've conducted a fairly extensive search and drawn mostly blanks. I found this http://confound.me.uk/maps/ppv4.pdf which looks like a good technique, but I'm unsure how it would work with multiple overlapping enclosed lines (areas).


texts: Why the icon can't display with my JSON data on Leaflet

texts:
I tried many time, but the icon didn't show on the map. I want to know what the problem on my code.

This is my code:

<!DOCTYPE html>

    <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.4/dist/leaflet.css"
          integrity="sha512-puBpdR0798OZvTTbP4A8Ix/l+A4dHDD0DGqYW6RQ+9jxkRFclaxxQb/SJAWZfWAkuyeQUytO7+7N4QKrDh+drA=="
          crossorigin="" />


texts: SpatialFilterClass to WFS filter for ArcObjects?

texts:
I want to translate a MapServer Query URL to a WFS GetFeature query.
Then I try to serialize a SpatialFilterClass to WFS filter.
any idea ?


texts:
I am working with a Digital Elevation Model for the area in and around a solar plant site.

The problem is that there are two DEM’S available with me, one with surveyed terrain (terrain inside of solar plant boundary with 1m resolution ) reflecting actual topography of site; simultaneously I have used SRTM DEM for topography outside of solar plant site(30 m resolution )

I think I could get better results of flood inundation mapping, if I am able to merge both elevation models (SRTM + Surveyed terrain ) by "re-instating" .

Does anyone have any experience doing this (or similar) and would like to share the method they used?

I believe it is possible to use the Raster Calculator to modify elevations in a DEM within a specific area using a mask in Q GIS, but I do not have experience with this kind of operation and I have not been able to find any example scripts or tutorials just from searching the web.


texts: Converting the result of pgr_dijkstra routing to an encoded polyline

texts:
I have this query for routing which is working, and would like it to give me an encoded polyline as a result instead of segmented geometries because I want to display it easily on a mobile Google MapView.
Select geom from pgr_dijkstra(
 'Select gid as id, source, target, st_length(geom) as cost from network', 1,3000, false) as di
 JOIN network pt
 ON di.edge = pt.gid ;

As I understand I must use ST_AsEncodedPolyline() function but I don't understand how to use the result of pgr_dijkstra() in a way it be used as parameter for ST_AsEncodedPolyline().


texts:
I'm am encountering an issue in ArcMap 10.2 regarding the Split Geoprocessing tool. I have a shapefile of roads (polylines) and another shapefile of districts (polygons). The roads shapefile is a result of running the MultipartToSinglepart tool so I can be absolutely sure there are no multi-part features.

I input the roads into the Split tool to be split on the districts. The problem is the output shapefile from the tool has not only split the lines based on the district boundaries, but has also divided my single-part roads into a multi-part feature. The divisions of these parts occur wherever the road intersected another road. As far as I can tell, multi-part features should not be created by the Split tool in this way. Does anyone know how to prevent the Split tool from creating multi-part features.

Tool is ran with these parameters, where m2s=roads:

An example of what I'm seeing is here:


texts: Split overlapping buffer polygons into zones

texts:
I am using ArcGIS 10.3 and
I need to calculate the mean altitude in some zones around specific points.
I used buffer tool for the desired distances from the point/center, but overlapping polygons have been created. I would like to keep only the zones between the predefined distances as individual polygons.
For example, I buffered at distances 2, 5, 10, 50, 150 meters and the corresponding polygons have been created. I would like to keep somehow individual polygons/zones/sections as following:
0-2m
2-5m
5-10
10-50
50-150
in order to calculate after the mean altitude in them.

Is there any way?

Preferably using toolbox procedures.


texts: Align line labels parallel to each segment of polygon offset from opposite ends

texts:
I have a polygon layer which I need to label on the opposite ends of a line and still maintain the labels parallel to this segment of a polygon. The labels need to be offset slightly from the end point of the line segment. In order to achieve this, I have derived lines that make up a polygon and now am struggling to achieve the desired labelling effect. I have used function $x_at( -1). How can I label it to be identical to

where I have the lines making up a polygon as


texts: opening vector layer in QGIS (3, 32bit) python (3.6) standalone application?

texts:
I know this question has been asked before, but the answer in Opening vector layer in QGIS Python standalone application? doesn't seem to work. i have tried both r"\path\to\file\" and "/path/to/file".

Im currently trying to open a shapefile/vector layer in a pyqgis standalone script. the script:

from qgis.core import *
import qgis.utils
import os

GUIEnabled = True
app = QgsApplication([], GUIEnabled)
app.setPrefixPath("C:/OSGeo4W/apps/qgis", True)
app.initQgis()

data_source = "C:/python_qgis/airport_tutorial/ne_10m_airports.shp"
layer = QgsVectorLayer(data_source, "ne_10m_airports.shp", "ogr")
if not layer.isValid():
    print ("Layer failed to load!")

It always return the "Layer failed to load!" statement. The shapefile has previously been opened in Qgis, and i have checked if the file is valid (it is).
Has something been changed in qgis version 3 that affected pyqgis?


texts: Shapefile appears correct on QGIS, but wrong on Google Earth

texts:
I've created a shapefile on QGIS, drawing over a raster that I've just georeferenced. The raster is a map with UTM lines drawn every 50 meters, so it was easy to georeference. Moving the mouse over it shows that the coordinates are accurate. The shapefile created also appears correct. (That's when the Project Properties' CRS is SAD69 / UTM zone 23S = EPSG:29193 -- for lat 20°S and lon 44°W -- so that I can see the coordinates in meters. The CRS of the shapefile drawn is WGS 84 = EPSG:4326, which automatically adjusts to the project, since 'on the fly' CRS transformation is allowed.)

When I open the same shapefile on Google Earth, I see it's shifted 50 meters to the West, when the difference between the Google Earth imagery and the map I used is just about 6 meters. What may be going on?

EDIT

The shapefile: https://www.sendspace.com/file/u22o3j

(Forgot to add a .prj file -- probably the source of all this). It contains:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

I did reproject (save vector layer as...) to SIRGAS 2000 (EPSG:31983), but got the same error. Tried to EPSG:29193, and now there's an error in another direction (north-nortwest, instead of west).


texts:
UPDATED
To clarify my previous post.  I have written a small application to convert maps saved as  PDF's  to PNG files then tile the maps to be used in a webservice.  After experimenting with the PNG Maps and updating my GDAL library to a newer version (1.9.2) Nuget GDAL in my program. I started experimenting with PDF's directly and not going through the time intensive task of converting them to PNG files.  A couple of things I have noticed in switching to PDF's is that:

1) the 256 pixel tiles are of better quality.

2) the tiling process runs much faster.

I have a lot of maps to convert on a daily basis, so the speed and quality improvements are a plus. One thing that pops up on quite a few PDF's that I tile is the exception: "Cannot Find GPTS Object".  This happens immediately when I go to open the dataset with the following code:

if(System.IO.Path.GetExtension(pathFile) == ".pdf" ||   System.IO.Path.GetExtension(pathFile) ==".png")
            {
                ds = Gdal.Open(pathFile, Access.GA_ReadOnly);// blows out on open
            }
            else
            {
                ds = Gdal.Open(pathFile, Access.GA_Update);
            }

I have searched quite a few hours on this exception and cannot find any documentation on it other than the source code that the exception is generated from.  I was wondering if anybody had run across this before and if there is a solution to the problem.

The PDF maps are not georeferenced.


texts: Altering layer's CRS to correct Layer Misalignment in QGIS

texts:
I've got a vector layer in my project which isn't aligned with the rest of the project (see attached JPEG).  These are work files so I'm not sure when they were created or how they might have been messed with since they were created.  We moved from ArcGIS to QGIS a couple of years ago.

The infilled blue vector layer is the one I need to move, slightly north and east, so that it correlates to the blue line (also a vector).  All layers (and the project) are currently in EPSG:27700 - OSGB 1936.

If I need to create a new CRS based on OSGB 1936 and give it some new parameters, then how would I go about this?

Further information:  This is on Exmoor in SW England, and the data is uniformly almost exactly 8m out to the SW point to point - i.e I need to shift it 8m NE.


texts: Joining z-values (from point file) to map generated by Kriging in ArcGIS Pro

texts:
I am using ArcGIS Pro 2.4.0.

I am making a fertilizer application map that changes the rates, based on the location of the field. The overall goal of the operation is to create a grid of squares, each with a grid containing a result of the Kriging interpolation.

I have a point file with nitrogen nutrient levels, the result of a soil-sampling operation:

The attribute table (with a row with for each point) that comes with it looks like this:

I want to know the specific level of nitrogen (NO3) in the field, so I used Kriging to interpolate it:

Now I have to convert from a floating-point raster to an interval raster. To do this, I first converted the floating-point raster (produced by Kriging) to an integer raster using the "Int" tool. Then, I used the "Raster to Polygon" tool to create my grid map. The grid produced is below:

I'm trying to get each of the grids to represent one of the Kriged nitrogen levels, but when I do a joins-and-relates action, I get a table with a bunch of unknown values.

I need the grids to reflect the values of the Kriging map that was generated earlier. How can I do that? Am I missing a joins-and-relates somewhere?


texts:
I built a server object extension in C# and Visual Studio 2010, using the ESRI sample for the Spatial Query SOE as a starting point.

(see the page on resources.esri.com)

The SOE never displays a property page in ArcCatalog or in Server Manager. If I enable it, it won't allow the map service to start, saying that the CLSID for the MapServer is not registered.

I think this is a bogus message, but I can't debug this - it doesn't even hit write any of my log messages to the server log. All the libraries I use (interfaces, extension, property page, server manager and register) have all been registered and GAC'd. The property page class shows up in Categories.exe and the assemblies all show up in the assembly folder.

It seems like there's some unstated step or steps that I haven't followed, but I can't find much on this. Has anybody actually built a working SOE? Any suggestions?


texts: Identifying spatially contiguous clusters in raster data using R

texts:
I am doing a classification of z over a single raster of dimension x,y,z.
r <- raster(nrows = 200, ncols = 400, xmn=-1, xmx=1, ymn=-1, ymx=1)
r[] <-  sample(1:100, size = 80000, replace = TRUE)

I want to do K-Means with Spatial Constraints:
Any of the solutions presented here would suit me
but my computer is too slow for that.
There must be a way to better leverage the raster format to make the process more efficient.
Any clue?


texts: Data stuck in arc seconds after being reprojected to UTM

texts:
I am working with world soil drainage data from the Harmonized World Soil Database v 1.2. The data started as an ascii file with no spatial reference and a resolution of 0.080595956, 0.080595956 Arc seconds, I later converted it to a raster file for ease of use.

After converting it to a raster I then defined the projection of the data into WGS_1984_World_Mercator as the data encompassed the world and I figured that's probably what the data was originally projected on (the original projection was not in the meta data).

The data projected fine then I proceeded clip the data to the study area and then reproject the data into WGS_1984_UTM_Zone_37S as this was the study area. The data reprojected but appeared on a different layer then the rest of my working data that is projected into WGS_1984_UTM_Zone_37S.

When I checked the source data the projection changed but the cell size appears to still be in Arc seconds:

I want the cell size units to be in meters. So far I have tried reprojecting and defining the data to make sure it wasn't an arc map error.

I also tried to reset the display units and the coordinates system from the data frame properties window. If you have any idea of what is causing this problem or have a solution let me know.


texts: Joining raster and vector layer using QGIS?

texts:
I have a raster layer (DTM of central Scotland) and a vector layer that shows the location of archaeological sites. I want to know how to join the layers so that I can figure out what height each archaeological site is at.

I have tried converting the raster layer into a vector layer and then doing a spatial joint but that does not seem to be working since when I polygonise the layer I seem to lose information.

Does anyone know how to "ask" QGIS what sites fall within what height or how to join a raster and a vector layer?


texts: Mapbox Android basics: import com.mapbox.mapboxandroiddemo.R

texts:
I am trying to rebuild the example code of the Mapbox first steps tutorial for Android Studio (https://docs.mapbox.com/help/tutorials/first-steps-android-sdk/). But I am struggeling with an import statement called

import com.mapbox.mapboxandroiddemo.R;

I get an error message in Andriod Studio: "Cannot resolve symbol 'mapboxandroiddemo'. I thought I would have to change that to my project name, but still getting the same error.

Whats wrong with this import statment and how to correct it?


texts: Using reader Feature Type Name as Writer Feature Type Name in FME

texts:
I've search and searched for a way to do this, but no method seems to work.  Basically I want to use the input feature class name (feature type in FME) as the output name.  This normally works, however the tool I'm using (ChangeDetector) adopts the "revised" input as the output Feature name instead of the "original" input name.  Even linking the "original" input name to the output name attribute does not work.  It keeps the "revised" input name.  I've tried using various transformers like the SubstringExtractor to pull the name and set the writer to fan out by this attribute, but even that fails.

While I would prefer to get this working within the workbench, I do plan on using this workbench in a python script, so it would be handy to get it working through that as well.  I setup the workbench to run using

subprocess.Popen(["pathToFme", "pathToWorkbench"])

This works just fine.  Subsequently I tried adding parameters at the end, to input the name based on code samples I found while searching this...

subprocess.Popen(["pathToFme", "pathToWorkbench"] --fileName: "nameOfOutput")

with "fileName" being the name of the parameter within the workbench, but I can't seem to get it to run, let along pick up the parameter input from the script.  I couldn't get os.system to properly execute the workbench, so I want to stay away from that.

Basically what I am trying to achieve is the ability to take a large number of inputs, run them through the changeDectector, and output with the same name as the "original" input.

Any help or suggestions would be appreciated.  Thanks.


texts: NameError: variable not defined even though set multiple times

texts:
This code original worked when run through the python command window in ArcCatalog. After I added the script to a toolbox I have had all sorts of troubles.

The most frustrating error has been:

NameError: name 'projectGDB' is not defined

Here is a section of the code. Unfortunately, the error persists. Again, the error only occurs when run as a tool from a toolbox. The NameError does not occur when copied and pasted into the command window.

# Import arcpy module & other helpful
import arcpy, os, xlwt

# Name Project
projectName = arcpy.GetParameterAsText(0)
if projectName == "#" or not projectName:
    projectName = "TEST"

# Create Project Folder
projPath = arcpy.GetParameterAsText(1)
if projPath == "#" or not projPath:
    projPath = "F:\\GIS\\TEST"
arcpy.AddMessage("Folder path for project is " + projPath)
print "Folder path for project is " + projPath

# Create Project Geodatabase
try:
    arcpy.CreateFileGDB_management(path, projectName, "CURRENT")
except Exception as e:
    print('Error creating gdb {0} at {1}'.format(projectName, projPath))
    print('Exception: {0}'.format(e.message))

# Set Work Environment
projectGDB = os.path.join(path, projectName + ".gdb")
arcpy.env.workspace = projectGDB
arcpy.env.scratchWorkspace = projectGDB
arcpy.env.overwriteOutput = True
arcpy.AddMessage("Workspace is set to " + projectGDB)

#Where is Landcover.gdb located?
landcoverGDB = arcpy.GetParameterAsText(2)
if path == "#" or not landcoverGDB:
    landcoverGDB = "F:\\GIS\\TEST\\Landcover.gdb"

# Load NRI Datasets
building = landcoverGDB + "\\imperv_building"

#[script continues for another 200 lines]

  Traceback (most recent call last):   File
  "F:\GIS\Packaged\NRI_sumAcres.py", line 23, in
      print "Project geodatabase is " + projectGDB NameError: name 'projectGDB' is not defined

Generally speaking, I do not understand why the script ran without error when copied and pasted into the python window of ArcCatalog, but now is dysfunctional once added to a toolbox.


texts: Change maxZoom option in runtime to ol.View in Openlayers 3

texts:
Im trying to change the maxZoom option once a functionality of the map is activated. So it has to be during runtime so it can be rolled back to the original maxZoom.

Creating the ol.View you configure this parameters like this:

new ol.View({
                ...
                zoom: 10,
                maxZoom: 17,
                minZoom: 10,
});

however, the api only allows to change zoom with setZoom()


texts: Finding Historical Maps of Soviet Union circa 1940?

texts: Anyone know where I can find a basic historical map of the Soviet Union, circa around 1940?