texts: Google earth image appears flat and in the wrong location when georeferencing

texts:
I've been trying to add an image of my site which I got from Google Earth Pro.  But I can't seem to georeference it properly using the georeferencing tool in Predictor noise modelling software.

I've followed the georeferencing tutorials and tried a bunch of things to try to sort it out but I've run out of ideas.  I've input the coordinates that I got from Google Earth,  and have gone through all the steps exactly but the image appears in the wrong location (next to my site), a lot bigger than it should, and completely flat. You can even see when grabbing the image and seeing the preview in the georeferenced tool that it isn’t correct.  It actually doesn’t grab the exact screen shot that I have ready with google earth.   I’ve also confirmed the coordinates in google earth pro match with the dxf that I’m using in the model.  For some reason the georeference tool imports it completely flat and in the wrong location, (wrong coordinates).

Have  you come across this before?


texts: OpenLayers ModifyFeature for multiple layers

texts: Is it possible to configure the ModifyFeature control to edit features in multiple layers?


texts: Why does georeferencing reduce the image resolution?

texts:
I'm trying to use a raster as a base image in my QGIS project.  It is a building schematic which is fairly high resolution.  Importing this into QGIS, it seems that the only way to re-position an image is through georeferencing.  If I do this, the resolution of the image drops significantly, to the point of now being usable for any presentation setting.

I'm trying to do one of two things, and I don't care which:

Reposition and rotate a raster layer; or
Georeference a raster layer without sacrificing image quality.


texts:
Basic question, but I have been asked to export data to an XYZ format. Is this the same as CSV files.

I am using QGIS 1.8.0. So i am presuming i select the data I want to export, right click its parent file in the window on the right and save selection as...CSV??


texts: Create shapefile from file.csv using python pyshp module

texts:
I followed the guide given here Automate the writing records from csv to shapefile using pyshp module. I follow this procedure and the code doesn't give me any error, but it doesn't produce any output. What can I do to get an output file.shp to load in QGIS?
My example file.csv is:

y,x,day1,day2,day3
37.0925854357,14.3860984446,0.0000079572,0.0000939491,0.0000028935
37.0978554286,14.5199850398,0.0000041594,0.0000821038,0.0000015805

This is the code I wrote:

import csv, shapefile
points = shapefile.Writer(shapefile.POINT)
with open("/path/to/file.csv", 'rb') as csvfile:
     csvreader = csv.DictReader(csvfile)
     header = csvreader.fieldnames
     [points.field(field) for field in header]
     for row in csvreader:
        points.point((float(row['y'])),float(row['x']))
        points.record(*tuple([row[f] for f in header]))
points.save('points')


texts: Is the 2010 Landsat world map (GLS 2010) available in a GoogleMaps browser?

texts:
Is the 2010 Landsat world map, also known as GLS 2010, available to view anywhere on the internet in a GoogleMaps browser?

I need to find it so that it can be added as a layer to a GoogleMaps interface on a wiki site Im helping out on.

It has much more up-to-date coverage in remote locations than Bing and Google, even though the resolution isnt so good.


texts: create temporary point geometry feature from coordinates in OGR with Python

texts:
I have two shapefiles containing two types of geometry features: one contains polygons and the other linestrings.

Because I'd like to use the Within method to check whether the end(or start) of a linestring falls within a polygon located in the other shapefile, I need to access the first and last points (as geometry features).
At the moment I iterate through the linestring feature and can identify the points using geometry.GetPoint(x) which returns the x and y coordinates of that point.

But when I'd like to use within to check whether a given point falls in the polygon object I can't because as Paolo explained in this question I need the geometry reference (using GetGeometryRef())of a point which can be retrieved like:

...
pointshp = driver.Open('test.shp')
pointlyr = pointshp.getLayer(0)
point = pointlyr.GetNextFeature()
point_geom = point.GetGeometryRef()
point_geom.Within(polygon_geom)

but I'm actually opening a shp that has linestring features. So after skimming through Utah State University GIS slides, I came up with this:

def within(pointcoord):
driver = ogr.GetDriverByName('test.shp')
ds = driver.CreateDataSource('test.shp')
layer = ds.CreateLayer('test',geom_type=ogr.wkbPoint)
fieldDefn = ogr.FieldDefn('id',ogr.OFTInteger)
layer.CreateFiled(fieldDefn)
featureDefn = layer.GetLayerDefn()
feature = ogr.feature(featureDefn)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(pointcoord[0],pointcoord[1])
feature.SetGeometry(point)
feature.SetField('id',1)
layer.CreateFeature(feature)
pointlyr = layer.GetLayer(0)
point = pointlyr.GetNextFeature()
point_geom = point.GetGeometryRef()
for x in xrange(len(polygons))
        point_geom.Within(Network.polygons[x].GetGeometryRef())

Which looks really inefficient to me. So just for the sake of using a method(within) I need to create not only a layer, but a shapefile, when all I really need is a simple point feature from which I can get the geometry reference to compare it with the one of the polygons.

The question is, how to create a point feature (temporary one, to be destroyed) from x and y coordinates? Or is there a simpler way of doing that than the example above?

On another hand, is there any other way of getting the geometry reference of the start or end point (as a geometry feature) of a linestring feature?


texts:
I'm trying to run a script that reads a file in a file geodatabase. I've been running the script in ArcMap 10.5 before without problem, so I think that might have something to do with it.

The script reads a file from a file database and uses in a Con and then save the Con layer as a file.

arcpy.env.workspace = r"C:\\Users\\...\\Output\\NinePoints.gdb"

rastList = arcpy.ListRasters()

for raster in rastList:
    print (raster)
    arcpy.env.workspace
    outCon = Con(raster, 1, "", "VALUE < 100")

    outputReclass = r'C:\Users\...\Output\Area.gdb' #lagres i ninePoints GDB'en
    outputFileName = "area"

    outRasterName = '{}{}'.format(raster, outputFileName)

    outCon.save(os.path.join(outputReclass, outRasterName))

The error I get is this one:

arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000865: Input conditional raster: flate2175reclasscost does not exist.
Failed to execute (Con).

It actually runs on the first file in the folder, but not on the second, which i find very weird.

I added the arcpy.env.workspace to the script and somehow it works now.


texts: Projection issues with Leaflet - polygon appears off

texts:
Using Leafletjs to upload coords from a json file as a polygon.
The original coordinates are Vicgrid94(3111).
I used QGIS to change the projection to LatLng (4326).

When I load the points onto my Leaflet map though, they're slightly off.
If I zoom in they get more accurate, but I want to limit the zoom.

My firewall at the office is blocking JSFiddle, sorry.

Here is my relevant code:

var map = L.map('map',{
  center: [-38.610980250630381,143.88599774132959],
  zoomControl: true,
  zoom: 13,
  maxZoom: 17,
  minZoom: 3
});

L.tileLayer('http://{s}.tile.openstreetmap.fr/hot/{z}/{x}/{y}.png', {
  maxZoom: 19,
  attribution: '&copy; OpenStreetMap, Tiles courtesy of Humanitarian OpenStreetMap Team'
}).addTo(map)

function loadData(count){
  thisShape = L.geoJson(data[count], {
    style: function (feature) {
      return {
        color: '#000',
        weight: 1
      };
    }
  }).addTo(map)
}
loadData(0)

This data is an array of geoJSON objects.

I tried placing crs: L.CRS.Simple and continuousWorld: true in the map and tilelayer objects, but this just changed the location of the map centre and polygon. Similar results for L.CRS.EPSG3857 and L.CRS.EPSG4326

EDIT

I've also tried exporting my data as projection 3857. But this did not seem to show up on my map, including when I set crs: L.CRS.EPSG3857

It seems the L.CRS.Simple projection wouldn't have worked anyway.

EDIT #2

Here is the same map, with a different geoJSON object, at a further level.
There's even more of the fire area in the ocean.


texts: Current field of view specific legends in QGIS

texts:
Is it possible to make QGIS generate a legend that is based only on the field visible in composer?

There was a ticket marked as resolved five months ago with this question but it doesn't divulge the solution!
http://hub.qgis.org/issues/3121

There are other questions appearing to ask this but I can't find the solution within them. A workaround I'm using at the moment is to save a new vector layer based on a selection of the field of view and plotting a legend from that.

Is there must be a more sensible way?


texts:
I set up GeoServer on one of my old laptops a few weeks ago with a connection to the internet over my home cable internet line (see: http://bit.ly/e9OUNP and http://bit.ly/hj1a6J)

I'm slowly learning the ins and outs of GeoServer, GeoExplorer, and other aspects of serving maps on a website. I assumed with this set up that I (or anyone else) would be able to access my GeoServer layers in ArcMap through WMS.

I installed the OpenGeoSuite Community Edition in the default directory. I have published layers with the proper settings and check boxes enabled. I can load the layers in a GeoExplorer map embedded in a web page (see: http://bit.ly/ecpGE2).

I have tried to load these layers into ArcMap using the following request (and many variations of it):
http://donmeltz.dyndns.org/geoserver/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetCapabilities&
but I cannot get this to work.

A couple of hours of searching has not turned up a solution to my problem yet. I'm hoping someone on here can point me in the right direction, and provide me with some ideas of what to look at next.


texts: Is the ESRI ArcGIS definition of California Albers projection using units of feet correct?

texts:
There is a map projection promulgated by the government of the state of California and recommended for use in certain state and municipal regulations.  It is an Albers Equal Area projection with origin at a convenient point amidst the state and appropriate false easting and false northing.  From what I can tell, the projection was first defined for use with linear units of meters.  It's listed in ArcGIS under State Systems, NAD 1983 California (Teale) Albers (Meters) with well-known identifier (WKID) 3310 from European Petroleum Survey Group (EPSG). At least, that's how it appears in the ArcGIS 10.4 catalog of projections.  The list of projections has other entries and WKIDs too for very similar projections.  Anyway, ArcGIS displays the list of parameters in the projection definition.  Among the projection parameters is a false northing value of –4000000 meters.  So far, so good.

The ArcGIS coordinate system list also includes one labeled NAD 1983 California (Teale) Albers (US Feet) with WKID 102600, on the authority of ESRI.  Units of feet are pretty common still in the US, so this projection sounds like a fine one I should use.  But when I examine the projection's definition and parameters, there's something odd.  The projection's linear unit parameter is the US foot.  No surprise there.  But the false easting value is –4000000.  That's –4000000 feet, I presume?

Shouldn't the false easting be in terms of US feet, which would be something like -13123333 feet?  (The factor is 3.28 or so, the meter-foot conversion).


texts: QGIS: "Wikimedia unlabelled layer" opacity

texts:
I'm using QGIS and I want to show a section of the basemap "Wikimedia unlabelled layer" (from the OpenLayers plugin) in full colouring.
But even if in the Properties box opacity is set to 100, still the map is very pale and colourless.
Is there a way to enhance colour quality? Maybe I have to use another similar OSM-file?


texts: Insert a Point into PostGIS using Python

texts: What is the right way to insert a Point to a PostGIS database using Python?


texts: Reading ellipsoid parameters using PyQGIS

texts: I would like to retrieve the semi-major axis of ellipsoid of a layer (any layer) using Python API. It seems that crs class can read only the ellipsoid acronym (?) Apparently, QgsDistanceArea class reads all ellipsoid parameters, but I don't know how to use it.


texts: How to prevent adding/editing with AttachmentEditor?

texts:
I have a WAB (for developers, version 1.4, JS API 3.16) widget with an AttachmentEditor. Depending on the user's privileges (managed through MVC), I want to prevent some of the users from modifying the attachments (this I can do by hiding the editor) and still showing them the attachments, but I found no way to let the editor show the attachments without permitting the user to add/delete them. Do I need to hack the editor code and use the tweaked version, or is there any clean option to do this? Bypassing AttachmentEditor completely seems to be the hardest way, so I hope it is not necessary.

I reposted this on Geonet.


texts: Convert NZTM to Latitude/Longitude in SQL

texts:
I have NZTM coordinate data stored in SQL Server in a Geometry field.

I'd like to convert these coordinates to Lat/Long so that it can be used with Bing maps in Excel (and a few other things).

I have found the formula to convert from NZTM to Lat/Long here on the LINZ website.

I'm not a mathematician, and that formula looks horrific.

Can anyone share a working example of the code in VBA or SQL that I could use to convert this on the fly?


texts: Filter nearby points leaving only one in QGIS

texts:
I have a shapefile with points, however the algorithm gave not one point per location, but several. I would need to filter these points, so that there is only one point, does not matter which one, from each group of points. The example of original points and points after filtration can bee seen below:

I know that this could be done using SAGA GIS Points filter option, however as I have 3.22 million of these points, the tool is not working. I also tried Dynamic point clusters in QGIS solution with QGIS Point displacement tool, however, at first it seems that points has been filtered leaving only one, but when zooming in, one point is really 3 points. Also this gives a wrong total point count.

Is there other option how this could be done using QGIS or other tools? Perhaps a python solution?


texts: Is there a way to properly have GDAL raise exceptions in Python?

texts:
I am working with GDAL in Python and have wanted to get GDAL to raise exceptions so that I can do something useful with them.  For instance if a file cannot be opened with gdal.Open() I would like to throw an exception instead of gdal returning None. I would just like a way to capture the errors instead of them being dumped to sys.stdout.

Currently I have tried using the gdal.UseExceptions() function but came across an issue.  When gdal.UseExceptions() is set it seems to hide certain errors and not raise an exception on them.  For instance ERROR 5: Access window out of range in RasterIO() Requested (1,15) of size 25x3 on raster of 26x17 disappears and no exception is raised. This seems pretty dangerous.

Any ideas are appreciated, thanks.


texts: Re-setting Radius of existing Circle object in esri/geometry/Circle

texts:
I would need to update the radius of existing circle object in  esri/geometry/Circle.

I've tried this code but doesn't work as expected.

var radiusSymbol = new SimpleFillSymbol(SimpleFillSymbol.STYLE_SOLID,
        new SimpleLineSymbol(SimpleLineSymbol.STYLE_DASHDOT, new Color([255, 223, 71]), 2)
        ,new Color([255,255,0,0.25])
      );

var pointX = new Point({
            longitude: X.longitude,
            latitude: X.latitude
        });
var circle = new Circle(
            pointX,
            {
                geodesic: true,
                radius: 2,
                radiusUnit: units.METERS
            });

var graphicRadius = new Graphic(circle, radiusSymbol);
map.graphics.add(graphicRadius);

Whenever new radius value and new lat&long are received, I need to update the circle with new radius value.

 var newCircle = new Circle(
         newPoint,
            {
            geodesic: true,
            radius: newRadius,
            radiusUnit: units.METERS
            });

    graphicRadius = new Graphic(newCircle, radiusSymbol);
    graphicRadius.setGeometry(newPoint);

It's not throwing errors but the circle is getting disappeared.

I hope there is setRadius method to use but it seems either I don't know how or not exist.

solved with graphicRadius.setGeometry(newCircle ,radiusSymbol);


texts: I would like to extract raster information (i.e. coordinate system, cell size) to excel file using python, how can I do that?


texts: Why is Openlayers KML-layer undefined when using permalink

texts:
I am working on Openlayers to dynamically load KML tracks on map. I want it to load filtered KML everytime map is moved or zoomed based on viewport and zoom level.

I got it working with code below, but for some reason Permalink doesn't work anymore. When pressing permalink link Debug console gives error: Uncaught TypeError: Cannot set property 'loaded' of undefined on line 108. That is the line layer.loaded = false; on UpdateKmlLayer() function.

Could someone help on what I am missing here?

// Start position for the map (hardcoded here for simplicity)
var lat = 65.37371
var lon = 24.49389
var zoom = 6;

var map; //complex object of type OpenLayers.Map
var KMLlayer;

//Initialise the 'map' object
function init() {

    map = new OpenLayers.Map("map", {
        controls: [],
        maxExtent: new OpenLayers.Bounds(-20037508.34, -20037508.34, 20037508.34, 20037508.34),
        maxResolution: 156543.0339,
        numZoomLevels: 16,
        units: 'm',
        projection: new OpenLayers.Projection("EPSG:900913"),
        displayProjection: new OpenLayers.Projection("EPSG:4326"),
        eventListeners: {
            "moveend": moveEnd,
        },

    });

    // Controls
    map.addControl(new OpenLayers.Control.Navigation());
    map.addControl(new OpenLayers.Control.PanZoomBar(), new OpenLayers.Pixel(0, 15));
    map.addControl(new OpenLayers.Control.Permalink());
    map.addControl(new OpenLayers.Control.ScaleLine({geodesic: true}));
    map.addControl(new OpenLayers.Control.MousePosition());
    map.addControl(new OpenLayers.Control.Attribution());

    // Osmarender
    layerTilesAtHome = new OpenLayers.Layer.OSM.Osmarender("Osmarender", {
        numZoomLevels: 18
    });
    map.addLayer(layerTilesAtHome);

    // Local KML
    KMLlayer = new OpenLayers.Layer.Vector("KMLlayer", {
        projection: map.displayProjection,
        strategies: [
            new OpenLayers.Strategy.BBOX(),
        ],
        protocol: new OpenLayers.Protocol.HTTP({
            url: "kml.php",
            format: new OpenLayers.Format.KML({
                extractStyles: true,
                extractAttributes: true
            })
        })
    });
    map.addLayer(KMLlayer);

    var switcherControl = new OpenLayers.Control.LayerSwitcher();
    map.addControl(switcherControl);
    switcherControl.maximizeControl();

    if (!map.getCenter()) {
        var lonLat = new OpenLayers.LonLat(lon, lat).transform(
            new OpenLayers.Projection("EPSG:4326"),
            map.getProjectionObject()
        );
        map.setCenter(lonLat, zoom);
    }
} //init()

function moveEnd(event) {
    console.log(event.type);
    UpdateKmlLayer(KMLlayer);
}

function UpdateKmlLayer(layer) {
    console.log("updateKML");
    console.log(layer);
    layer.loaded = false;           // setting loaded to false unloads the layer
    layer.setVisibility(true);      // setting visibility to true forces a reload of the layer
    layer.refresh({ force: true, params: { 'zoom': map.getZoom()} });   // force new KML data
}


texts: Land use change detection with large amount of vector files

texts: Im doing a project where I'm creating a map that shows how land use changed from 1990-2000 for the city of Minneapolis. I have three files, one of the cities boundaries, one with the land use of the 2000's and the other has the land use of 1987, 1990, 1997, and 2000. I merged and intersected the land use .shp with the city land cover .shp. However I'm having trouble showing how land use changed. I'm having difficulty because the land use history .shps have over 80,000 files. Because of the large number of files its not easy for the naked eye to see the polygons. Any ideas on how I can make the changes more easily perceivable.


texts: Out_Folder_Path as a Parameter and OS.Path.Join

texts:
First, I am new to StackExchange and to Python.

I would like to allow my user to select the location to store results of a TableToTable conversion.

I do not understand how to convert the out_folder_path and out_name to parameters. My understanding is that the out_name is part of the out_folder_path. If that is the case, how can out_folder_path exist before out_name is created?

Can I set the out_folder_path equal to the workspace? (I have successfully set the workspace as a parameter)

I have read several posts regarding using os.path.join.

Below is my attempt to add t he code following the directions in advanced help.

#import arcpy
import arcpy

#import OS
import os

#setup environment
from arcpy import env

#set overwrite to true
env.overwriteOutput = True

#Set workspace as parameter
arcpy.env.workspace = arcpy.GetParameterAsText (0) # -Type = workspace - let's user select database of shapefiles to convert to tables/spreadsheets

#Variables for creating folder
##Syntax: CreateFolder_management (out_folder_path, out_name)
out_folder_path = arcpy.GetParameterAsText(1) # Type = workspace; allows user to select location to save results
out_name = arcpy.GetParameterAsText(2) #Type = string; allows user to name the folder for storing results

#Create A Folder to Store Spreadsheets
arcpy.CreateFolder_management(out_folder_path, out_name)

#Get a list of shapefiles in the folder
fcs = arcpy.ListFeatureClasses()
fcs_count = len(fcs)

for fc in fcs:
  outfc = os.path.join(out_folder_path, fc)
  ##Syntax: TableToExcel_conversion (Input_Table, Output_Excel_File, {Use_field_alias_as_column_header}, {Use_domain_and_subtype_description})
  arcpy.TableToExcel_conversion(fc,outfc)


texts:
Trying to find the best way to take the current map view extent

const extent = mapView.calculateExtent();

and get a WKT from that value. i thought i might be able to take that extent, convert it into a polygon:

const poly =  fromExtent(mapExt);

and then somehow get the geometry from that polygon and convert to wkt like so:

const wktFormat = new WKT();
const polyGeom = //get geometry from poly here?
const ret = wktFormat.writeGeometry(polyGeom, {
  dataProjection: 'EPSG:4326',
  featureProjection: 'EPSG:3857'
});

but i'm not able to find a good way to do that.  Am i missing a step?


texts: Error on Clip Raster By Mask Layer (QGIS)

texts:
Can not seem to figure out this problem when i am running Clip Raster by Mask Layer function:

>
INFO    gdalwarp -of GTiff -cutline "C:/Users/user/Downloads/gis files/final 23.05.2020.shp" -cl "final 23.05.2020" -crop_to_cutline crs=EPSG:3857&format&type=xyz&url=https://tile.openstreetmap.org/%7Bz%7D/%7Bx%7D/%7By%7D.png&zmax=19&zmin=0 "C:/Users/user/Downloads/gis files/tst.tif"
2020-05-23T18:09:16     INFO    GDAL execution console output

         FAILURE: No target filename specified.

         Usage: gdalwarp [--help-general] [--formats]

          [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid]

          [-order n | -tps | -rpc | -geoloc] [-et err_threshold]

          [-refine_gcps tolerance [minimum_gcps]]

          [-te xmin ymin xmax ymax] [-tr xres yres] [-tap] [-ts width height]

          [-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]

          [-srcnodata "value [value...]"] [-dstnodata "value [value...]"] -dstalpha

          [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]

          [-cutline datasource] [-cl layer] [-cwhere expression]

          [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]

          [-of format] [-co "NAME=VALUE"]* [-overwrite]

          [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*

          [-doo NAME=VALUE]*

          srcfile* dstfile

         Available resampling methods:

          near (default), bilinear, cubic, cubicspline, lanczos, average, mode, max, min, med, Q1, Q3.

         Required parameter missing -

         The system cannot find the file specified.

         'url' is not recognized as an internal or external command,

         operable program or batch file.

         'zmax' is not recognized as an internal or external command,

         operable program or batch file.

         'zmin' is not recognized as an internal or external command,

         operable program or batch file.


texts:
I want to have a loop on a date variable in Google Earth Engine, but it hangs and show me nothing when running just below simple code:

var aDate = ee.Date.fromYMD(2015,12,21);
var bDate = ee.Date.fromYMD(2015,12,31);
for(var d = aDate; d < bDate; d.advance(1,'day')){
  print(d);
}

I tried it in different forms and also via while loop, but it always hangs.
Any clue?


texts: Join attributes by location gives out wrong sum

texts:
I have an input vector polygon layer which has different areas.
Now I want to join a different vector polygon layer which contains over 1 million 100x100 meter squares with information about the population.

As a result I want to see which part of the squares belongs to which area of the other polygon and also have the sum of the population of each area.

Everything is working correctly, except the sum-function: The layer with the squares has a combined population of little more than 80 million. But after the join the result has either 76 million or 84 million depending which kind of geometric predicate im using. But never the correct 80 million.

Both have the same CRS and I'm using QGIS 3.4.3 Madeira.


texts:
I would like to identify the Kth order neighbors of an edge on a network, specifically the neighbors of a large set of streets.  For example, I have a street that I'm interested in looking at, call this the focal street.  For each focal street I want to find the streets that share an intersection, these are the first order neighbors.  Then for each of those streets that share an intersection with the focal street I would like to find their neighbors (these would be the second order neighbors), and so on...

Calculating the first order neighbors using ArcGIS' geoprocessing library (arcpy) took 6+ hours, second order neighbors are taking 18+ hours.  I have posted on this in the ESRI forums HERE (link includes arcpy code that works but is crazy slow).  Needless to say I want to find a more efficient solution.  I have created a python dictionary which is keyed on each street and contains the connected streets as values.  For example:

st2neighs = {street1: [street2, street3, street5], street2: [street1, street4], ...}.

street 1 is connected to street 2, 3, 5; street 2 is connected to street 1 and 4; etc..  There are around 30,000 streets in the study area most have fewer than 7 connected streets.  A pickled version of the data used in the code below IS HERE.

I assumed that knowing the first order neighbors would allow me to efficiently trace the higher order neighbors.  But the following code is providing incorrect results:

##Select K-order neighbors from a set of sampled streets.
##saves in dictionary format such that
##the key is the sampled street and the neighboring streets are the values

##################
##IMPORT LIBRARIES
##################

import random as random
import pickle

#######################
##LOAD PICKLED DATA
#######################

seg_file = open("seg2st.pkl", "rb")
st_file = open("st2neighs.pkl", "rb")
seg2st = pickle.load(seg_file)
st2neigh = pickle.load(st_file)

##################
##DEF FUNCTIONS
##################

##Takes in a dict of segments (key) and their streets (values).
##returns the desired number of sampled streets per segment
##returns a dict keyed segment containing tlids.
def selectSample(seg2st, nbirths):
    randSt = {}
    for segK in seg2st.iterkeys():
        ranSamp = [int(random.choice(seg2st[segK])) for i in xrange(nbirths)]
        randSt[segK] = []
        for aSamp in ranSamp:
                randSt[segK].append(aSamp)

    return randSt

##Takes in a list of all streets (keys) and their first order neighbors (values)
##Takes in a list of sampled  streets
##returns a dict of all sampled streets and their neighbors.
##Higher order selections should be possible with findMoreNeighbors
##logic is the same but replacing sample (input) with output from
##findFirstNeighbors

def findFirstNeighbors(st2neigh, sample):
    compSts = {}
    for samp in sample.iterkeys():
        for rSt in sample[samp]:
            if rSt not in compSts:
                compSts[rSt] = []
            for compSt in st2neigh[rSt]:
                compSts[rSt].append(compSt)

    return compSts

def findMoreNeighbors(st2neigh, compSts):
    for aSt in compSts:
        for st in compSts[aSt]:
            for nSt in st2neigh[st]:
                if nSt not in compSts[aSt]:
                    compSts[aSt].append(nSt)
    moreNeighs = compSts
    return moreNeighs

#####################
##The nHoods
#####################

samp = selectSample(seg2st, 1)
n1 = findFirstNeighbors(st2neigh, samp)
n2 = findMoreNeighbors(st2neigh, n1)
n3 = findMoreNeighbors(st2neigh, n2)

#####################
##CHECK RESULTS
#####################
def checkResults(neighList):
    cntr = {}
    for c in neighList.iterkeys():
        cntr[c] = 0
        for a in neighList[c]:
            cntr[c] += 1
    return cntr

##There is an error no streets **should** have 2000+ order neighbors
c1 = checkResults(n1)
c2 = checkResults(n2)
c3 = checkResults(n3)

Help!


texts: Perfect circle buffer on QGIS Madeira 3.4

texts: Working on a road junction I would like to make a perfect rounded circle buffer of 1 km around this junction in QGIS 3.4. Here is what is look like now.


texts: How to increase or decrease the spatial resolution of a band in Google Earth Engine?

texts:
In LANDSAT

var mndwi = ee.Image("LANDSAT/LC08/C01/T1_RT_TOA/LC08_148044_20130425").normalizedDifference(['B3', 'B6']).rename('MNDWI')

gives us MNDWI Band

In Sentinel 2 the resolutions of Green Band (B3 is 10m) and SWIR Band (B11 is 20m) are different.
How to change SWIR Band to be finer resolution so that I can apply normalized difference?


texts: Styling multiple features in shapefile using GeoServer?

texts:
I have a shapefile loaded in a GeoServer instance, and I need some guidance on how to best style the output. The file is depicting public/private land ownership in the State of Colorado. The shapefile has different feature properties (adm_name) for each agency - State Land, BLM, US Forest Service, Private, etc. Each type of ownership needs to be a different color.

A PNG image of the unstyled map is here.

I need it to look like this (without the grid, of course).

So, this is as much a theoretical question as a technical one.

Can I just apply the styling with a SLD file? Or inline in the web browser through the Leaflet code? Or should the features be split up into different layers and a single style applied to each layer? Or something else completely?

It's currently a huge resource to load on a web page. So, any suggestions to accomplish this efficiently are most welcome.


texts: ArcMap Intersect Tool is missing just one intersect point of two polyline shapefiles?

texts:
I'm trying to find the intersect of 2 polyline shapefiles. One shapefile is a long line (of multiparts), the other shapefile is 11 straight lines that "fan out" radially. The ends of the first shapefile are snapped to around the midpoint of the first and last of the 11 lines, (i.e. the two radially-farthest lines) and crosses each of the other lines around their midway.

I need to know the exact X & Y location of the intersect points. I've been running the intersect tool, then "add XY". There should always be 11 intersect points (because one of the shapefiles is 11 straight lines), but it is repeatedly missing the "last" point (radially farthest to the right, or clockwise) and this messes up some data processing later on down the line (even if I manually add the point to the intersect points shapefile).

I've tried:
- running "repair geometry" first - I don't actually know what this does, but it didn't help
- converting both shapefiles to "singlepart" (because as is, the first consists of multiple parts of each line segment that I mapped, and the second is 11 lines) then running the intersect on the singlepart shapefiles
(neither worked)

does anyone have any idea why this doesn't work, and what I can do to fix this?


texts: Using QGIS Processing-Toolbox in stand alone Python scipt (Microsoft Visual Studio)

texts:
I'm using QGIS 3.14, Python 3.7, Processing-Version 2.12.99 and Microsoft Visual Studio 2019 (IDE).
Everything works fine, the only problem I have is the import of the processing-algorithms. Related to the answer of @Germán Carrillo: Using QGIS3 Processing algorithms from standalone PyQGIS scripts (outside of GUI) I have:
import sys

from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)
from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath(r'C:\OSGeo4W64\apps\qgis', True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append(r'C:\OSGeo4W64\apps\qgis\python\plugins')

import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

When it comes to line "import processing" I get the following error "triggered exception: no module named '_uuid'":

I tried different paths/variables etc. but nothing works. I did not come across a similar problem during my research.


texts: How to create an if else condition with ArcGIS Javascript API for slider and imageLayers from switch statement

texts:
The code below is the switch statement that I am currently using.  I can't seem to get the values (case) to match up with the horizontal slider I am using (4 unique values-each layer).  I was thinking of trying an if else condition instead.  Anyone have the best way to accomplish this?

function setPixelFilter() {
        var val = slider.get("value");
        switch (val) {
          case 0:
            var imageLayer = new ArcGISImageServiceLayer("https://...ImageServer")
            map.addLayer(imageLayer)
            break;
          case 100/3:
            var imageLayer = new ArcGISImageServiceLayer("https://...ImageServer")
            map.addLayer(imageLayer)
            break;
          case 100/3*2:
            var imageLayer = new ArcGISImageServiceLayer("https://...ImageServer")
            map.addLayer(imageLayer)
            break;
          case 100:
            var imageLayer = new ArcGISImageServiceLayer("https://...ImageServer")
            map.addLayer(imageLayer)
            break;


texts: What is the most efficient way to get a list of featureclasses from a large ArcSDE database?

texts:
I have a SDE DB which consists of 40 feature datasets and around 300 featureclasses.  I need to get the full list of all featureclasses whose name does not contain "ARCHIVE". It takes several minutes using the script below. It will be unacceptable if I publish this as a geoprocessing service. How can this be improved?

import arcpy,os
import json
from arcpy import env

env.workspace=r"C:\ysde\lcimrl.sde"
arcpy.AddMessage("--------get resource dataset \n")

resourceFeatureList=[]
lodFeatureList=[]
dsList = arcpy.ListDatasets("*","Feature")
for ds in dsList:
    dsDesc = arcpy.Describe(ds)
    dsName=dsDesc.name.upper()
    arcpy.AddMessage("--------get feature class list from "+dsName+" \n")

    fcList = arcpy.ListFeatureClasses("*","",ds)
    print fcList
    for fc in fcList:
        desc = arcpy.Describe(fc)
        fcName=desc.name.upper()
        #print dsName+"_"+fcName
        if   "ARCHIVE" not in fcName:
            if dsName!="SDE.ALTERNATIVES":
                resourceFeatureList.append(fc)
            else:
                lodFeatureList.append(fc)
print resourceFeatureList
print lodFeatureList


texts: Retain "region" elements for when using fortify() on a SpatialLinesDataFrame

texts:
I'm layering multiple plots using ggplot2 to create a map of Alaska that combines school location and the major road network. I've learned that despite being labeled "major roads" there is wide variance in quality, which is actually encoded in the original shapefile data in a variable called MIN. I would like to illustrate road quality by coloring the line segments, but I'm unsure how to retain data elements from the original SpatialLinesDataFrame object after fortify() in order to use them later to change the color or shape of the lines.

This is a simplified version of the map I would like to make with just two layers of the plot (the state outline and roads):

To replicate my exact circumstance you can download my data here.  The road data has a variable called MIN that seems to encode road quality, which I would love to join to the final fortified data frame:

library(ggplot2)
library(rgeos)
library(maptools)
library(rgdal)

coast <- readOGR(".", layer = "Alaska_Coast_1000000_py")
proj4string(coast) <- CRS("+init=epsg:3338")
coast <- spTransform(coast, CRS("+proj=longlat +datum=NAD83"))
coast.df <- fortify(coast, region = "TYPE")
# fix Aleutian Islands
coast.df$long <- ifelse(coast.df$long>0, coast.df$long*-1, coast.df$long)

roads <- readOGR(".", layer = "Major_Roads_ln")
proj4string(roads) <- CRS("+init=epsg:3338")
roads <- spTransform(roads, CRS("+proj=longlat +datum=NAD83"))

# where I would love to have a "region" option as with polygons:
roads.df <- fortify(roads)

# this plot doesn't try to color the lines,
# but is the basis for what I'm hoping to adjust
a4 <- ggplot() +
      geom_polygon(data = coast.df,
                   aes(x = long, y = lat, group = group),
                   color = "gray", size = 0.25, fill = "white") +
      geom_path(data = roads.df, aes(x=long, y=lat, group=group),
                size = 1)
# could imagine adding color = "road_type"

print(a4)


texts: Is there a way to download more than 13MB of OSM data to QGIS?

texts: I'm using QGIS to make maps with OSM streetmaps data. Problem is that with QGIS -> Vector -> OSM -> download data, after about 13MB I get "Download error Operation cancelled". Sometimes I get "Download successful" after 8MB but can't import the .osm because of "XML error, Premature end of document". I reduced the areas of the download to keep the size down and two out of five were incomplete files. With Maperitive I've downloaded bigger files, from rambler.ru.  Is there some setting in QGIS to change the source of the data? Or some other way of solving the problem?


texts: Sentinel 2 surface reflectance error Python code Google Earth Engine

texts:
We are trying to do the surface reflectance correction using Python codes in Spyder. The study area is defined as a polygon geometry and it's a fairly large area. We want to get all sentinel images that covered the study area in 2017 and correct all of them.

However when we run the python code, error popped up saying there is too many pixels in our study area:

EEException: Error in map(ID=20170103T172712_20170103T173023_T13QGE):
Image.reduceRegion: Too many pixels in the region. Found 220713615, but only 10000000 allowed.

We tested - the pixel number 220713615 should be the pixels of the geometry we specified (the study area). We tried a smaller geometry with less than 10000000 pixels, and the code worked perfectly.

Weird things are: we are not using any "Image.reduceRegion" function in our code, and it doesn't seem to have anything to do with the pixels of the study area if we want to surface correct the satellite images.


texts: Creating and Labeling time series (features)

texts:
I am trying to carry out crop type identification using optical remotely sensed data of Sentinel-2A. I am new to both fields (remote sensing and machine learning) and I want to use time series (NDVI, NDWI...) as my input data (training), and I also have a vector shapefile as my labels.
How can I create these data time series?
I am able to create an NDVI map (per raster image/date), but I was not successful with creating a time series per pixel using SNAP. (this is a pixel-based classification)


texts: ArcMap raster calculator raising negative values to a power

texts:
I have a raster calculator expression:

Power("dataset.tif",1.5)

where dataset is float, and has some negative values.

When I run this expression in raster calculator, all of the negative values within the dataset have an output of 'nodata'.

How can I overcome this issue?


texts:
I wanted to know how could I do a offset line of one layer in a new layer.

I take my layers, even if I only have one active (streets)

canvas = iface.mapCanvas()
layers = canvas.layers()

In my case, layers is just one layer, and this layer have 600 lines.
After that I don't know first create a new layer where put the result (the offsets) and the other problem create the offset ( -2 and 2). I saw other post trying to resolv my problem, but I tried to do that but didn't work.

I did that manually, but I want that in python.

Manually, I had this:

And after offset I got this:

But I tried to write several codes and all codes didn't work.


texts: Disable antialiasing in GeoServer .png output?

texts:
I would like to completely disable antialiasing of polygon borders in Geoserver's WMS png output. The reason is because we do downstream colour remapping, and intermediate colours break that.

This is the problem:

I would like that to look sharp and jaggy :)

EDIT

I see from the documentation that there is a "vendor parameter" that should allow this:

&FORMAT_OPTIONS=antialiasing%3Aoff

Alas, it does not work for me. (Geoserver 2.6.4)

EDIT 2:

Here's a live URL you can try (for a while...)

Changing the DPI changes the fidelity of the rendered geometry, but doesn't affect antialiasing.


texts: MODIS emissivity extraction using ArcGIS Desktop?

texts:
I  am doing validation of 40 Landsat images.now I need two spectral emissivities of the 8-day MODIS land surface temperature and emissivity product (MOD11A2).and I have to extract emissivities for MODIS band 31 and 32 for each imaging date.

How do I do it with ArcGIS?


texts:
I installed GSL 1.16 & GDAL 1.10, then I proceed to install QGIS 1.8.0

After running the .pkg(installer) I get the following message.

Failure
The GDAL 1.9 framework is required.

Does this mean I should be using GDAL 1.9 instead of 1.10 or is it just not recognizing 1.10?


texts: How to create a separate environment using ArcGIS Pro Python in Visual Studio

texts:
I am switching from Pycharm to Visual Studio as my Python script IDE.

I have been trying to create an environment for ArcGIS Pro but I had no success so far.

I can create Custom Environment from the existing environments but it directly uses the environment for ArcGIS Pro.

I would rather use the ArcGIS Pro Python (3.6) as a base interpreter and create a separate environment.

Anyone has the same problem or am I doing something wrong here?


texts:
I have 3 collared animals and their GPS position is mostly daily recorded, included information about the date, time, longitude, latitude. The daily records vary from any record for weeks to 24 records per day. I manually reduced records for 3 records per day as much standardized as possible. I would like to know how to calculate the daily average distance they have traveled even though the data are very unbalanced.
I am using QGIS 2.18.12.


texts: Convert shapefile to CSV including attributes AND geometry?

texts:
I have a shapefile with 60k+ entries, all of which are polygons with corresponding attributes (acreage totals, landowner names, tax ID #s, etc.).  What I ultimately need is a CSV file with all of these attributes and their corresponding geometry (in the KML compatible xyz format, that is, NOT the WKT format).

I know that I can open the .dbf file in Excel and get the attributes.  I also know that I can open the shapefile in QGIS and copy the data into Excel, which gets me attributes and WKT geometry.

Is there a simple way to convert the shapefile to CSV (openable in Excel) with attribute and Google Earth friendly geometry?


texts:
I am trying to make a model using ArcGIS Model Builder which would iterate the select (analysis) tool.

Basically, I have a feature class with more than 10,000 points. These point have different Z values (height), with the field name is "POINT_Z".
So, I want to export these points into different feature class by dividing them based on the Z values.

Without using model builder, I will need to manually select these points using the select(analysis) tool.
This is the SQL script to select the features that have height between 32 and 33

(POINT_Z >=32) AND (POINT_Z <33)

This process of using select(analysis) tool is repeated with different height range, but I was wondering whether there is a way to iterate this process.
The lowest and highest values of the POINT_Z is 32.355 and 69.804 respectively.

Thus, the intended output is 38 different feature classes altogether (FC with Z value between 32-33, FC with Z value between 33-34, so on and so forth)
Can anybody help me on this problem?


texts: Why doesn't animated show up as an option in Carto Builder Style Tab?

texts:
I'm trying to do a time series in carto.

My data table has the following:

-polygon geometry  (geom)
-polygon name (string)   - eg State name
-category (string)
-date  (datetime)

Someone asked the question I'm looking to answer as Is it possible to add a timeline play button in CARTO Builder?, but even though I have a timedate column, and have been able to add the widget, when I look under the style tab, I don't have the animated option.

Does anyone have any idea what the issue is?


texts: Autocomplete "select" not centering on map

texts:
Everything seems to be working good until I select from the list of names.  The map isn't centering on the proper coordinates.  I'm using JQuery UI and and OpenLayers 3.  here's my code:

$(function () {
    function log(message) {
        $("
").text(message).prependTo("#log");
        $("#log").scrollTop(0);
    }
    var centerCoords;
    $("#city").autocomplete({
        source: function (request, response) {
            $.ajax({
                url: "http://ws.geonames.org/searchJSON?username=XXXX",
                dataType: "jsonp",
                data: {
                    featureClass: "P",
                    style: "full",
                    maxRows: 12,
                    name_startsWith: request.term
                },
                beforeSend: function () {
                    $('.input-group-addon').html("");
                },
                success: function (data) {
                    response($.map(data.geonames, function (item) {
                        return {
                            label: item.name + (item.adminName1 ? ", " + item.adminName1 : "") + ", " + item.countryCode,
                            value: item.name + (item.adminName1 ? ", " + item.adminName1 : "") + ", " + item.countryCode,
                            lat: item.lat,
                            lng: item.lng
                        };
                    }));
                    $('.input-group-addon').html("");
                }
            });
        },
        minLength: 2,
        select: function (event, ui) {
            log(ui.item ?
                    "Selected: " + ui.item.label :
                    "Nothing selected, input was " + this.value);
            var location = ol.proj.transform([ui.item.lng, ui.item.lat], 'EPSG:4326', 'EPSG:3857');
            console.log(location);

            var a = map.getView();
            a.setCenter(location);
            a.setZoom(14);

            var marker = new ol.Feature({
                geometry: new ol.geom.Point(location)
            });
            var vectorSource3 = new ol.source.Vector({
                features: [marker]
            });
            vectorSource3.addFeature(marker);

            var vectorlayer3 = new ol.layer.Vector({
                source: vectorSource3
            });

            var iconStyle = new ol.style.Style({
                image: new ol.style.Icon({
                    anchor: [0.5, 46],
                    opacity: 0.75,
                    scale: 0.10,
                    anchorXUnits: 'fraction',
                    anchorYUnits: 'pixels',
                    src: "../images/flag.png"
                })
            });
            marker.setStyle(iconStyle);
            map.addLayer(vectorlayer3);

        },
        open: function () {
            $(this).removeClass("ui-corner-all").addClass("ui-corner-top");
            $(".ui-autocomplete").css("z-index", 1000);
        },
        close: function () {
            $(this).removeClass("ui-corner-top").addClass("ui-corner-all");
        }

    });
});


texts: Would it make sense to translate a city to latitude=0 longitude=0 before doing any calculations?

texts:
I'm working on a conversion system for Unity. We want to import some geo-referenced data (not models) losing the least floating point resolution possible.
I am currently working on the idea of moving all the data as we import it right to the equator across the Greenwich line, in other words to the geodetic origin.
Basically we take the commonly accepted center of that town and translate it.
This would simplify the conversion to cartesian space and assure that we can stay within acceptable bounds for floating point resolution (usually ±10000).
Does this all make sense?


texts: ArcObjects IConflictsWindow3 not seeing version reconciliation conflicts

texts:
We are working on a reconcile and post add-in for one of our departments. I'm able to successfully complete the task if there are no conflicts or if I specify which version 'wins'. We would like to utilize the integrated ArcMap functionality for conflict display and resolution.

Unfortunately, when trying to display the conflict window, I'm having issues. First of all, the conflict window does not specify that there are any conflict classes even though the IVersionEdit4 object specifies there are and the reconcile result is 'true'. Secondly, if I try to display the conflict window anyway, the window opens and ArcMap crashes--even if I surround the code in a try block.

I'm open to suggestions. I am not particularly interested in manually finding the sources of the conflict when the functionality is there to accomplish the task in ArcMap natively and the code is exposed in ArcObjects to accomplish the task.

Am I missing something? Do I need to somehow assign which class the conflict is to the window?

Here's the code in question:

    bool ReconcileToDefault(Version ToReconcile, out bool? hadConflicts)
    {
        ChangeDataSourceToDefault(ToReconcile); //changes datasources in map document to default version to prevent lock issues
        hadConflicts = null;
        IVersionedWorkspace wkspc = (IVersionedWorkspace)DBConnect.ConnectToDB("sde.DEFAULT", true);
        if (wkspc == null) return false;
        IVersion Reconcile = wkspc.FindVersion(ToReconcile.name);
        //EventListener evtListener = new EventListener(Reconcile);
        IMultiuserWorkspaceEdit muWorkspaceEdit = (IMultiuserWorkspaceEdit)Reconcile;
        IWorkspaceEdit workspaceEdit = (IWorkspaceEdit)Reconcile;
        IVersionEdit4 versionEdit = (IVersionEdit4)muWorkspaceEdit;
        muWorkspaceEdit.StartMultiuserEditing(esriMultiuserEditSessionMode.esriMESMVersioned);
        workspaceEdit.StartEditOperation();
        //Open the conflict resolution window. (ArcMap bombs out here, so I'm trying a different way)
        //if (cWin == null)
        //{
        //    IEditor editor = (IEditor)ArcMap.Application.FindExtensionByCLSID(new UID() { Value = "esriEditor.Editor" });
        //    cWin = (IConflictsWindow3)editor.FindExtension(new UID() { Value = "esriEditor.ConflictsWindow" });
        //    cWin.Reset();
        //    cWin.Visible = true;
        //}
        hadConflicts = versionEdit.Reconcile4("sde.DEFAULT", false, false, false, true); //I set acquire locks to true if I'm actually posting.
        if (hadConflicts == true)
        {
            versionEdit.ConflictClasses.Reset();
            var cl = ConflictResolutionDialog.CurrentClass; //variable is null with current database state
            IConflictClass first = versionEdit.ConflictClasses.Next(); //variable is a valid IConflictClass with current database state
            //MessageBox.Show(ConflictResolutionDialog.HasConflicts().ToString()); //displays false
            //MessageBox.Show(String.Format("Had {0} Conflicts.",ConflictResolutionDialog.ClassCount)); //conflict count = 0
            //ConflictResolutionDialog.Visible = true; //causes arcmap to crash
        }
        //versionEdit.Post("sde.DEFAULT"); //works just fine
        workspaceEdit.StopEditOperation();
        workspaceEdit.StopEditing(false);//generally set this to true if I post
        return true;
    }

Thanks for your attention,


texts: Drawing Sectors / Conic Sections in OpenLayers?

texts:
I am having trouble finding working examples in OpenLayers 5.3 where conic sections and/or sectors can be drawn.

Does anyone have any working code or can point to the correct classes to use?

Here's some examples of shapes I am looking to draw

https://www.mathworks.com/help/examples/map/win64/CreateSouthernHemisphereConicProjectionExample_01.png
https://www.inchcalculator.com/wp-content/plugins/uber-calc/assets/img/shapes/sector.png

I am able to make a basic version of the shape by doing this (though I'm sure there's a way more efficient way):

const outer = fromCircle(new Circle(map.getView().getCenter(), myRadius), 128);
const half = Math.ceil(outer.getCoordinates()[0].length / 2);
outer.setCoordinates([outer.getCoordinates()[0].splice(0, half)]);

const inner = fromCircle(new Circle(map.getView().getCenter(), myRadius / 10), 128),
outer.appendLinearRing(new LinearRing(inner.getCoordinates()[0]));

const myFeature = new Feature(outer);


texts:
I'm looking to optimize point proximity geo searches time.

My input is lat,lng point and I'm searching on a precomputed set of locations to n nearest points.

I don't care how much time/space the building of the precomputed index of locations will take but I do care the queries will be super fast.

I'm thinking about using geohash as the search key, where I would first check if I get results for X chars of the key and then continue to trim down chars from the end of the key until I start to see results.

To my (very sparse for now) understanding of geo index techniques this approach should be able to produce the fastest results (in terms of query time) compared to all other known implementations (such as R Tree and co.)


texts: TopologyException side location conflict reading MapInfo with Sharpmap

texts:
I'm using SharpMap 0.9 to read various MapInfo (.tab) files.

I have used a number of files without any problem and been able to perform intersection queries on them like so:

var interesectionResults = provider.ExecuteIntersectionQuery(geometry);

where the provider is created like so:

var layerLocation = new FileInfo(_layerDirectory.FullName + "\\" + layerName +FileExtension);

if (!layerLocation.Exists)
    throw new Exception(string.Format("Could not locate the layer data file, {0}", layerLocation));

return new Ogr(layerLocation.FullName);

The geometry used in ExecuteIntersectionQuery is defined as follows:

var cadastreProvider = GetLayerProvider(Settings.Default.PropertyLayerName);

//Query cadastre to get property bounding box
var query = String.Format("SELECT * FROM {0} WHERE {1} = {2}", cadastreProvider.LayerName, Settings.Default.PropertyKeyField, propertyID);

var propertyResults = cadastreProvider.ExecuteQuery(query);
var resultsTable = propertyResults.Tables.First();
var propertyFeatureRow =  resultsTable.Rows[0] as FeatureDataRow
var geometry = propertyFeatureRow.Geometry

For a particular set of new MapInfo files I have recevied for a client everytime I do the intersection query the following error is thrown

System.ApplicationException was caught
  Message=TopologyException: side location conflict 402348 6.9506e+006 1.7e-308
  Source=ogr_csharp
  StackTrace:
       at OSGeo.OGR.Layer.GetNextFeature()
       at SharpMap.Data.Providers.Ogr.ExecuteIntersectionQuery(FeatureDataSet ds) in E:\Investigation\GIS\Sharpmap\Sharpmap src\Trunk\SharpMap.Extensions\Data\Providers\OgrProvider.cs:line 483
       at SharpMap.Data.Providers.Ogr.ExecuteIntersectionQuery(Geometry geom, FeatureDataSet ds) in E:\Investigation\GIS\Sharpmap\Sharpmap src\Trunk\SharpMap.Extensions\Data\Providers\OgrProvider.cs:line 465
       at SharpMap.Data.Providers.Ogr.ExecuteIntersectionQuery(Geometry geom) in E:\Investigation\GIS\Sharpmap\Sharpmap src\Trunk\SharpMap.Extensions\Data\Providers\OgrProvider.cs:line 735
       at IFMA.Common.GIS.MapInfo.Implementation.GetLandAttributesForLayerByPropertyID(String propertyID, String layerName, Boolean displayLayerDetails, String layerIdentifierField, String layerDescriptionField, Int32 queryType) in E:\Common\Current\01- Trunk\GIS\Implementations\IFMA.Common.GIS.MapInfo\Implementation.cs:line 145
       at IFMA.MasterPlan.BusinessLogic.PropertySystemManager.GetApplicationPropertiesLandAttributes(Application application, Nullable`1 systemVersion) in E:\MP\Current\03- Working\IFMA.MasterPlan.BusinessLogic\PropertySystemManager.cs:line 211
  InnerException:

I can't seem to find much information about this error but I gather its at least being thrown by the underlying provider (OGR) that Sharpmap uses to read MapInfo files.

Has anyone got an idea what might be causing this?


texts: ArcMap mxd has strange artifact that won't delete?

texts:
I have an MXD with all our data organized and displayed with the appropriate symbology, that took quite some time to make. It is used as the basic map that a number of people use here and customize for their own purposes. However, there is a strange "artifact" that seems to exist in one spot. It is an old polygon of some data that has been since updated. It will show up in the map regardless of whether all the data is turned off or not. Even if all the data is removed, it still persists. I tried selecting it, thinking maybe it was some label turned to annotation but the "delete" in the right click menu is greyed out. It seems deleting the data frame will remove it, but remaking this mxd would be more trouble than just dealing with the artifact.

How can I get rid of this?


texts: How to properly Drape vectors over a DEM?

texts:
I have a DEM of my study area, and the streams in it.
Whenever I add both to ArcScene, and drape the vectors over the DEM, some of the vectors disappear from over the surface, and are visible under the surface. (see Image below).

Are there any specific settings to properly drape the vectors over the DEM?
Are there any other options or software that produce a better result (visually) when you have to drape vectors over a DEM?


texts: pgrouting - Choosing an algorithm with a predominance for a specific type of road

texts:
I would like to use one of pgrouting's algorithm to find the (shortest) path between 2 points with predominance for a specific type of road.

To build the data I used the osm2pgrouting.

Is there an algorithm where I could define my preferences for the road type (in table ways class_id column class_id which represents road type or in table osm_way_types)?  If not how should I tweak my database? Knowing that I used the osm2pgrouting --> there is table public.osm_way_classes there is max_speed but it seems poor choice for tweaking it. I would like to leave it alone. I didn't really find anything useful in other tables. Any ideas?

@edit (second question deleted) I saw that in dijkstra algorithm there is cost in inside query. Maybe add column priority to way_classes (or way_types). And then join the table ways with table way_classes and then with table way_types. So in query there won't be ST_Length / length in cost but length * priority. What do you think?


texts: ESRI JavaScript 4 - Shouldn't the editor widget honor definition expression?

texts:
If a feature layer has a definition expression defined, as in the code pen example, shouldn't you only be able to see and edit (select for editing) features that meet the definition expression?

It seems that the widget tool doesn't honor the definition expression and its quite confusing.

I need a work around.

Example: You have a feature layer in a web map, there are two features in the data source that over lap, the definition expression is true for 1 and false for the other (so you can't see the second in the web map),  if you select the edit tool and click on the one that is visible, both are queried from the feature service and both are displayed in the edit widget results.

Is there something I'm not configuring correctly?

Here is a code pen example.  Click on the pink zigzag patrol line with the edit tool active

https://codepen.io/AliWieck/pen/ExaMWJq?editors=1000

Line 91 of this codepen example applies a definition query that should limit display of any police route lines like so

layer.definitionExpression = "PatrolType = 'Segway patrol'"

yet when I click on a location where I know there are Segway and other types of routes that are overlapping, I get all features returned  and displayed in the editor widget even though only one  is visible on the map.  Click the button for any feature that doesn't meet the definition expression does nothing.


texts: What's wrong with my expression to control label colors?

texts:
I am struggling to alter the label colours according to field values (as strings) in the attributes table.

In my mind what I want is:

If the string in column X = 'Proposed'
Then label the point using the string in column Y, with color A

Else
Label the point using the string in column Y, with color B
End

I have tried a bit, but struggling with the syntax

case when  "ownerenc" = 'Proposed'
then  (color_rgb(256,256,256))
else  color_rgb(56,56,0)
end

But I realise this is quite wrong! Many thanks for your suggestions.


texts:
I am new to Geonetwork but my understanding is that after harvesting the Catalog has a set of (CSW) Records that remain 'flat' and there s no hierarchy.

In particular, i have a number of Products that belong to different Catalogs and i need to query the CSW endpoint and get back a number of Products that belong to Catalog 1.

How can i set their association?


texts: subquery in select by attribute in Arcgis

texts:
I have a table with several fields. field A contains the same kind of values as does field B. I tried to select all points where the value in field A does not exist in anywhere field B using query builder in select by attribute.

For example

ID A  B
1  a  b
2  b  c
3  c  a
4  d  b

It should select the points with ID 4 since value d does not exist in field B.
I think I should use a subquery, but how?


texts:
I'd just like to add some GIS layers, such as shapefiles in the Maps.me app.

Has someone already done it?


texts: Counting total number of other link ends link is connected to in QGIS?

texts:
I have a polyline layer (pretty much an OSM network) that I need to count the total number of other link ends the primary link is attached to. In other words, I have over 59,000 links in a city network and I'm trying to measure connectivity. I want to count the number of endpoints to other links that each link is attached to.

I also need to maintain the edge IDs, which rules out the sDNA package.

I have a feeling there are python scripts I could use, but I'm not familiar with python (intermediate with R, but unfortunately no exp with igraph).

It would be great if there was a qgis version of Lines   and Junctions tool from ArcGIS.


texts: Defining area for QuickOSM query in QGIS Modeler

texts:
I would like to use the QGIS Modeler to collect all objects of as OSM feature within polygons in a layer. The QuickOSM algorithm has a function called "build query inside an area" and I was hoping I could select my polygon layer here but that doesn't work (see picture below, "Inside the area").
Does anyone have an idea how to use this algorithm or if it is possible at all to select within polygons?


texts: Change color of feature (active cursor row) with ArcPy

texts:
I have a selection of polygon feature and with a cursor loop I walk through them.
I determine the extent of eachfeature and arcmap will zoom to the feature. My default selection color in ArcMap is light-blue. So the whole selection is light-blue.
After determining the extent of each feature I create a png-file on disc.
Is there a way that the corresponding feature for the selected cursor row will be red before I create my file on disc.


texts:
All of the Geopackages stored on one of my drives (D:/) appear to be corrupted, and don't work, but work fine if they are copied to another drive.   The other file formats on the same drive all work fine (shape files, tiffs, csv, etc). Any geopackage on my drive D:/ that I've tried to open with QGIS cannot be opened by other software (eg R) until it has been moved to another location.

The error also crops up when I try creating a new geopackage in the drive, but can create them in other locations. QGIS creates an empty geopackage, but then says it's not a valid data source and will not add it to the map.

I can create a geopackage using R on my drive D:/ and it works fine until I try opening it with QGIS, at which point it is not openable by either program.

I updated QGIS from version 3.2.0 to 3.2.3 after this problem showed up, but it did not change the behavior.

Does anyone know why Geopackages would stop working in one location, and how to fix it?  Using a different location isn't a great solution because the drive D:/ is the primary drive I use for GIS projects.


texts:
How can I draw a flying/moving line on Cesium?

I need a base line and a little point flying from one origin coordinates to a destination, and it needs to be animated.

Cesium only provide basic apis like drawing a line on the surface or towards a certain xyz point.


texts:
How to correct the extent of a raster?

According to the metadata of the raster I downloaded, it should be (WGS84):

North : 85.044
South : -85.044
East : 180
West : -180

But this is the extend of the raster when I display it in ArcGIS:


texts: Checking if ArcGIS extension has already been checked out manually in ArcPy?

texts:
I've written a python script that uses a Spatial Analyst function and hence requires the Spatial Analyst extension to be enabled.  I would like to be able to check out the extension only if it hasn't already been manually checked out by the user (i.e. via the "Customize->Extensions" dialog).

There doesn't seem to be any way of testing whether the extension is already checked out out using arcpy.CheckExtension as that returns Available regardless of whether a license is already checked out or not (assuming there is a license available).

My code can check out the extension and then check it back in easily enough which is fine if the code is running outside of ArcGIS. However, this causes an issue in ArcGIS if a user has already checked the extension out manually via the "Customize->Extensions" dialog as they will suddenly start getting the "ERROR 000824: The tool is not licensed" or "ERROR 010096:  There is no Spatial Analyst license currently available or enabled." errors when trying to use other Spatial Analyst tools/functions (as my code has checked the extension back in).

A further complication that may confuse the user is that when they open the  "Customize->Extensions" dialog to check if the SA extension is checked out, the tick will still be in  the extensions check box as checking the extension in via code does not update the dialog. E.g see screenshot below

My work around was to test using a simple spatial analyst tool (Int(1)) in a try:except: clause (see below), but that's just ugly.

from arcpy.sa import Int

def sa_is_licensed():
    try:
        Int(1)
        return True
    except (ExecuteError, RuntimeError) as e:
        # "ERROR 000824: The tool is not licensed"
        # "ERROR 010096:  There is no Spatial Analyst license currently available or enabled."
        if 'ERROR 000824' in e.message or 'ERROR 010096' in e.message:
            if arcpy.CheckExtension('Spatial') == 'Available' and arcpy.CheckOutExtension('Spatial')=='CheckedOut':
               return False
            else:raise
        else:raise

already_licensed = sa_is_licensed()

# My big long script...

if not already_licensed: arcpy.CheckInExtension('Spatial')

Are there any better ways to check from ArcPy if an ArcGIS extension has already been checked out?


texts: Point to line interpolation using heading in R

texts:
I asked this question at StackOverFlow, but got no answers, so I thought I would try here also:

I can draw a spline through 10 successive points like done in this post

set.seed(1)
n <- 10
x <- runif(n)
y <- runif(n)
p <- cbind(x,y)
xlim <- c(min(x) - 0.1*diff(range(x)), c(max(x) + 0.1*diff(range(x))))
ylim <- c(min(y) - 0.1*diff(range(y)), c(max(y) + 0.1*diff(range(y))))
plot(p, xlim=xlim, ylim=ylim)
text(p, labels=seq(n), pos=3)
xspline(x, y, shape = c(0,rep(-1, 10-2),0), border="red")

But is it possible to force the direction of the line in each point to be at a specific heading? - For example, if the heading was 0 in point #2, it would be forced to go straight up, and not as now in a north-easterly direction.

Like on a compass: 0=North, 90=East, 180=South and 270=West.

For example you could use the following 10 headings:

h <- round(runif(n, min = 0, max = 359),0)

Also, the output should be able to be saved as a spatial line.


texts: Debugging NameError: name 'QgsLinearlyInterpolateDiagramRenderer' is not defined from PyQGIS?

texts:
I have gone through an online tutorial on how to do a pie chart map using PYQGIS.

The following is my code, identical to what was provided.

# -*- coding: utf-8 -*-
from PyQt4.QtGui import*
from PyQt4.QtCore import*
layer = QgsVectorLayer("C:\Users\Aaron\Documents\pygis\County10PopnHou\County10PopnHou.shp","Population","ogr")
categories = [u'PCT_WHT',u'PCT_BLK',u'PCT_AMIND',u'PCT_ASIAN',u'PCT_HAW',u'PCT_ORA',u'PCT_MR',u'PCT_HISP']
colors = ['#3727fa','#01daae','#f849a6','#268605','#6810ff','#453990','#630f2f','#07dd45']
qcolors = map(QColor,colors)
diagram = QgsPieDiagram()
ds = QgsDiagramSettings()
ds.font = QFont("Helvectica", 12)
ds.transparency = 0
ds.categoryColors = qcolors
ds.categoryAttributes = categories
ds.size = QSizeF(100.0,100.0)
ds.sizeType = 0
ds.labelPlacementMethod = 1
ds.scaleByArea = True
ds.minimumSize = 0
ds.BackgroundColor = QColor(255,255,255,0)
ds.penWidth = 0
dr = QgsLinearlyInterpolateDiagramRenderer()
dr.setLowerValue(0.0)
dr.setLowerSize(QSizeF(0.0,0.0))
dr.setUpperValue(2000000)
dr.setUpperSize(QSizeF(40,40))
dr.setClassificationAttribute(6)
dr.setDiagram(diagram)
layer.setDiagramRenderer(dr)
dls = QgsDiagramLayerSettings()
dls.dist = 0
dls.priority = 0
dls.xPosColumn = -1
dls.yPosColumn = -1
dls.placement = 0
layer.setDiagramLayerSettings(dls)
label = QgsPalLayerSettings()
label.readFromLayer(layer)
label.enabled = True
label.writeToLayer(layer)
if hasattr(layer,"setCacheImage"):
    layer.setCacheImage(None)
layer.triggerRepaint()
QgsMapLayerRegistry.instance().addMapLayer(layer)

Upon running on QGIS using the Python console, I get this error.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/Users/Aaron/Documents/pygis/UsingPieChartSymbols.py", line 21, in
    dr = QgsLinearlyInterpolateDiagramRenderer()
NameError: name 'QgsLinearlyInterpolateDiagramRenderer' is not defined

I have no idea why this is occuring and what is missing.


texts: How to exclude some layers from ArcReader legend?

texts: I have created an ArcReader map for staff which has a number of imagery layers, such as satellite imagery and a raster topographic map. I also have a legend which automatically updates based on which layers are selected in the table of contents. How do I exclude the imagery layers, though, because they are not useful in a legend? When they appear in a legend they are just RGB colours.


texts: How to segmentize/sample a LineString using geodjango (postgis/geos) and/or openlayers?

texts:
I'm trying to sample points along a linestring. The goal is to have the distance between 2 consecutive points to be at most X (eg. 50m).

I'm drawing linestring in openlayers and sending it to a geodjango app. I can do the "segmentization"/sampling in the client (openlayers) or in the server (geodjango: DB API or Geos API). I don't know yet where I will do it, so I'm interested in both solution (if they exist).

I've searched but have not found any simple solution to this really simple problem (at least, it looks very simple).

I've seen this one for example: How to interpolate GPS Positions in PostGIS but it uses complex postgis query to do the interpolation.

Isn't there any 'split(nseg)' function somewhere that would split a segment into nseg sub segments ? Something close to what shapely offers ? (http://gispython.org/shapely/docs/1.2/manual.html#linear-referencing-methods ). I could use shapely, but that would add another dependency (not really justified...).

EDIT: I guess I should have a look at linearref from GEOS (https://code.djangoproject.com/ticket/11948 for the geodjango part)


texts: ArcGIS calculate difference between rows in the same field that is date time format

texts:
I'm looking for help tweaking the python code discussed here: Attribute Table Subtract Value From Previous Value

I'm having trouble figuring out how to adjust it for a date field in which I want to subtract a time by the previous time to get gaps.  Here is what I have so far and the error received:

import arcpy

from arcpy import da

vertices ="CBR05_T2_10_1_2014_WLanduse_1"
gap1=0
with arcpy.da.UpdateCursor(vertices, ["INDEX_", "Time", "GapDiff"]) as cursor:

    for row in cursor:
        gap2 = row[1]
        row[2] = datetime(gap2 - gap1)
        gap1 = gap2
        cursor.updateRow(row)

Runtime error
Traceback (most recent call last):
  File "", line 8, in
TypeError: unsupported operand type(s) for -: 'datetime.datetime' and 'int'


texts:
QGIS 2.18.11

This is not a question as such, but certainly something which I would learn from if resolved. I have uploaded a shapefile of polgyons from multiple layers all merged into one. I have made sure the features are singleparts, checked validity of the geometry (which checks out) but still spits out an empty output. I've been trying to select half the features and dissolve twice to reduce the pressure put on the algorithm. I have tried v.dissolve but the same result occurs when using vector>geoprocessing tools>dissolve...just an empty output.


texts: Model Builder buffer output not being created?

texts: Created a model that creates a buffer and uses that buffer to select by location features within the buffer and run various analysis.  My issue is...the buffer output ('QQ 1500 2000ft Buffer Output' below) is not actually being created as a feature class in the folder I specified.  The model still runs the analysis successfully, just doesn't create the actual feature class.  Only way to create it was to add the 'Copy Features' function.  Furthermore, I tried connecting the copied features output ('QQ Buffer Output' below) to processes and the feature again wasn't created.  So it seems to be an issue with connecting processes to the output that stops the create feature from working.  Anyone have any insight into this issue?  To note, the output feature class is based on model parameters so my folder path looks like '%ProjectFolder%\%FeatureName%.shp'


texts: Points Layer to Distance Matrix QGIS 3.6 - Noosa

texts:
I'm trying to calculate a sample distance matrix from four points in the BNG CRS:

Using Vector-> Analysis Tools -> Distance Matrix, I get:
Input point layer is a MultiPoint layer - first convert to single points before using this algorithm

How can I get a distance matrix from my layer of many points?


texts: Extracting and transposing a single record in db manager

texts:
This seems like a simple operation so apologies if I've missed an existing answer.

What I need to do is extract single record (a client) from a database of clients (in Spatialite) and present the query results as a transposed table.

So, to simplify, the client database would have fields; name, address, phone number, etc. for a large number of clients (one row each). I want a query which will return a table as below for a particular client I pick.

field,value

name, {value in 'name' field for selected client}

address, {value in 'address' field for selected client}

phone number, {value in 'phone number' field for selected client}

etc.

Hope that makes sense.


texts: Appending multiple FeatureClasses to one another gives ERROR 000735?

texts:
I have 5 fairly large Feature Classes with identical fields but different attribute data which I have been trying to merge/append together unsuccessfully.

When using merge_management the process runs successfully, but I am left with a feature class with many new fields created for each of the 5 feature classes (e.g. name1, name12, name123, name1234).

I really need to just retain the name field and add to it rather than creating new fields for the imported Series').

When I attempt to append rather than merge the operation fails with an

  ExecuteError: Failed to execute. Parameters are not valid. ERROR
  000735: Input Datasets: Value is required Failed to execute (Append).

Find my code below:

(For Merging attempt)

arcpy.env.workspace = r"D:/some.gdb"

fc_list = arcpy.ListFeatureClasses()

fc_list = [x for x in fc_list if "Pymade" in x]

print("Feature classes to join: " + str(fc_list))
print("Joining borough FCs into final output")
arcpy.Merge_management(fc_list, "MapPlutoAllBoro_{0}".format(today))

(For Appending attempt)

arcpy.env.workspace = r"D:/some.gdb"

arcpy.CreateFeatureclass_management(r"D:/some.gdb", "MapPlutoAllBoro_{0}".format(today))

fc_list = arcpy.ListFeatureClasses()

fc_list = [x for x in fc_list if "Pymade" in x]

print("Feature classes to join: " + str(fc_list))
print("Joining borough FCs into final output")
arcpy.Append_management(fc_list, "MapPlutoAllBoro_{0}".format(today))


texts:
I am used to working with MapInfo when using TAB files, however am in the process of learning how to use QGIS instead. One snag I have hit is that when rows in a table are deleted, the empty rows still exist within the table (as can be seen easily by opening the TAB in MapInfo). You cannot see the empty rows in QGIS but they are there.
Usually, in MapInfo, one would simply go to Table > Maintenance > Pack Table...
My question is: can you pack a table in QGIS?
I am using QGIS Desktop 3.10.13
[I have found that saving the file as a new file gets rid of the empty/deleted rows, but this is not practical]


texts: QgsVectorLayer.IsValid() always returns False?

texts:
I'm writing a standalone Python script which involves some geoprocessing using qgis libs. Unfortunetely, I couldn't manage to properly instantiate a vector layer - every time i run .isValid() on QgsVectorLayer instance, it returns False. I tried using both Shapefile and PostGIS layers.
I'm aware of the following threads:
* https://stackoverflow.com/questions/41449252/python-qgis-layer-failed-to-load
* Problem loading vector layer (pyQGIS)
but none of given answers resolves my problem.

Here is the output from QgsApplication.showSettings() :

Application path not initialized
Application path not initialized
Application path not initialized
Application state:
QGIS_PREFIX_PATH env var:               C:/OSGEO4~1/apps/qgis-ltr

Prefix:
Plugin Path:
Package Data Path:      C:/OSGEO4~1/apps/qgis-ltr
/.
Active Theme Name:
Active Theme Path:      C:/OSGEO4~1/apps/qgis-ltr
/./resources/themes\\icons/
Default Theme Path:     :/images/themes/default/
SVG Search Paths:
User DB Path:   C:/OSGEO4~1/apps/qgis-ltr
/./resources/qgis.db
Auth DB Path:   qgis-auth.db

And here is my code snippet

from qgis.core import *
import os
import sys

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = 'C:\\OSGeo4W64\\apps\\Qt5\\plugins'
os.environ['PATH'] += 'C:\\OSGeo4W64\\apps\\qgis-ltr\\bin;C:\\OSGeo4W64\\apps\\Qt5\\bin'

QgsApplication.setPrefixPath("C:\\OSGEO4~1\\apps\\qgis-ltr", True)
qgs = QgsApplication([], False)
qgs.initQgis()

sys.path.append('C:\\OSGeo4W64\\apps\\qgis-ltr\\python\\plugins')
import processing

project = QgsProject.instance()

# Gminy.shp is a valid shapefile (opens in QGIS), containing a .shx file as well
tshp_layer = QgsVectorLayer("G:\\PRG_7.05.2019\\Gminy.shp", "test_lyr", "ogr")
# this returns false
print(str(tshp_layer.isValid()))

# connection to PostgreSQL DB using proper credentials works fine
uri = QgsDataSourceUri()
uri.setConnection("<host>", "5432", "<db_name>", "<user>", "<pass>")
uri.setDataSource("<schema>", "<table>", "<geometry_field>")
drogi_layer = QgsVectorLayer(uri.uri(False), "my_custom_layer_name", "postgres")
# this also returns false
print(str(drogi_layer.isValid()))

qgs.exitQgis()

print('All done!')

Which returns:

False
False
All done!


texts: postgis: cluster/group points based on distance

texts: We have a set of thousands of points and want to group all of those which are within 100m distance from each other (to get the centroid from each group). A first idea was to build 100m buffers around each point and to union/group all intersecting/overlapping polygons to one (but there can also be non-overlapping polygons within one cluster, see screenshot) - after that, we could simply calculate the centroid from the merged polygons. But I do not know how to do that in Postgis (in QGIS there is a simple tool called 'dissolve' that merges overlapping polygons, but we want to automatize the process)


texts: Dividing value of polygon in polygon inside it

texts:
I'm working with population data. I want to divide the population value (labeled in the orange square)in the polygons inside it (blue polygons).

Any ideas?


texts: Accessing Attribute table in ArcMap on Mac computer?

texts:
I'm using ArcGIS for Desktop on a Mac computer via Parallels.

How do I access the Attribute Table?

I cannot right-click on the layer
I keep getting Layer properties...


texts: Transform coordinates in csv file from wgs84 to nad27 using qgis

texts: I've got a .csv file which lists boreholes and their x & y coordinate information in NAD27. I was wondering if there is a way to use QGIS when transforming the coordinate information into WGS84. The desired end product would be a .csv file where all boreholes are specified in the wgs84 system.


texts: Unexpected error merging three point shapefiles

texts:
I get the following error when I try to merge the files:

The error says that there is an unexpected error. I have deleted the lock files and an empty ArcMap document. It still gives me this error.


texts:
I'm experimenting with this new ArcGIS Pro 2.0+ that became available at my workplace (currently using ArcMap), and I'm experiencing some inconsistencies with datasets brought in via query layers from our SQL Server (currently '08 server. Tried w/ a '16 server and had same issues).

Some datasets seem to work just fine (selecting / "identifying" / searching), and some datasets erroneously return incorrect records when I "identify" them with the new "Explore" button. Example: I "identify" a utility pole in front of our building, which I know has a name of "POLE-X", but the "identify/explore" window returns back "POLE-W". If I "identify" that pole again it may return another pole, "POLE-Z", or sometimes no object attributes at all.

Has anybody encountered odd behavior like this before?


texts: Interpolate (TIN) within a polygon extent QGIS/GRASS

texts:
I am trying to interpolate (TIN) in QGIS to produce surface models. I need to interpolate within a polygon area extent. The options available in QGIS do not allow this as far as I am aware. Is there any way of doing this, including utilising GRASS, and if so, how do you do it?

The image below shows an example of a desired extent with the points lying within it.


texts: using GDAL installed using homebrew from python on OSX

texts:
I would like to work with jpeg2000 files in python but when I install GDAL using pip or conda it can't read jpeg2000s.

GDAL's version is: 2.1.0
module 'osgeo.gdal' from '/Applications/anaconda/lib/python3.5/site-packages/osgeo/gdal.py'

I found this guide to installing gdal 2 on osx which I followed and I now have a version of gdal installed which can deal with jpeg2000.

user$ gdalinfo --formats | grep -i jp
  JPEG -raster- (rwv): JPEG JFIF
  JP2OpenJPEG -raster,vector- (rwv): JPEG-2000 driver based on OpenJPEG library
  JPEG2000 -raster,vector- (rwv): JPEG-2000 part 1 (ISO/IEC 15444-1), based on Jasper library

My question is how can I make python use this version of GDAL and not the one that is downloaded by pip or conda with the gdal package?


texts:
I am trying to show a particular kind of rendering to my geometry when a user clicks on the record of a feature data grid and zooms to the selection, but can't figure out what I am doing wrong...

Here is my code for the featureLayerFromMap.UpdateCompleted += new System.EventHandler(fslayer_UpdateCompleted);

void fslayer_UpdateCompleted(object sender, System.EventArgs e)
    {
       ESRI.ArcGIS.Client.FeatureLayer featureLayerFromMap = (ESRI.ArcGIS.Client.FeatureLayer)(_mapService.Map.Layers["FeatureLayerControlLayer"]);

        //*** Flashing of Polyline and Polygon when zoomed from a feature data grid record
        _objectIdDic = new Dictionary<string, string>();
        int count = 0;
        featureLayerFromMap.Graphics.ToList().ForEach(x => _objectIdDic.Add(x.Attributes["OBJECTID"].ToString(), count++.ToString()));

        if (featureLayerFromMap.Graphics.Count > 0)
        {
            if (featureLayerFromMap.Graphics[0].Geometry.GetType() == typeof(ESRI.ArcGIS.Client.Geometry.Polygon))
                featureLayerFromMap.Renderer = this._mapLayout.Resources["PolygonGraphicLyrSelectedRenderer"] as IRenderer;
            else if (featureLayerFromMap.Graphics[0].Geometry.GetType() == typeof(ESRI.ArcGIS.Client.Geometry.Polyline))
                featureLayerFromMap.Renderer = this._mapLayout.Resources["PolylineGraphicLyrSelectedRenderer"] as IRenderer;
        }
    }

The feature layer is provided from a selection in a combo box... Can someone guide me on this?


texts: How to export MapServer and featureserver links from QGIS

texts: QGIS has an option to export wmf, wcs, wfs ad xyz links all in one go. I have a lot of MapServer and featureserver layers, how do I export those all in one go?


texts: OGR layer intersection TopologyException: Input geom 1 is invalid

texts:
I'm traying to intersec two shp files, for most cases it works, but some times it fails.

ERROR 1: TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 270434.99999640981 4298084.9998855386 at 270434.99999640981 4298084.9998855386

The code is more or less for two existing shp poDS1 and poDS2

ly1=poDS1->GetLayer(0);
ly2=poDS2->GetLayer(0);
ly3= poDS->CreateLayer(intersectFile,ly1->GetSpatialRef());
err = ly1->Intersection(ly2,ly3, NULL,ProgressFunc, (void *) "Intersection");

The result of  err is 6;
I can intersect it in qgis without err or advise.

Files are here where
76_marina.shp is poDS1 and
76_marina_LC82020332016182LGN00_R5 is poDS2


texts: Loop stops before expected without error message

texts:
I try to run a script that clips all the feature class of a geodatabase by the successive features of a clip layer (named pochoir).

After Accessing geometry with SearchCursor in ArcPy and copy only selected with CopyFeature_management, the test of the intersection is know working.

In the following code, all the "for row in rows" is working. The only issue is that the script stops just after processing the first fc in fcList.
For the first fc processed, there was three rounds. The two first were True for the intersecting test while the third was not. The two first were well clipped and the third was well ignored. But I think something block after the if or the else statement because next fc is not run.

I used try/except but no error were recorded.

Here is the code :

arcpy.env.workspace = r"Z:\Documents\SIG\tests\synthese_tocorrect.gdb"
fcList = arcpy.ListFeatureClasses()
pochoir = r"Z:\Documents\SIG\tests\decoupe.shp"#Clip layer with multiple entities
rows = arcpy.SearchCursor(pochoir)
outFolder = r"Z:\Documents\SIG\tests\sorties"
count = 0
for fc in fcList:
    for row in rows:
        secteur = row.secteur
        where_clause = "secteur = '{0}'".format(secteur)
        arcpy.MakeFeatureLayer_management(pochoir, "ptmp", where_clause) #Create current clip layer
        arcpy.MakeFeatureLayer_management(fc, "toselect")#Create tmp layer to use in the select by location (no feature class allowed)
        arcpy.SelectLayerByLocation_management("toselect", "INTERSECT", "ptmp")# testing the intersect
        outName = "{0}_{1}_{2}".format(fc, secteur, count)

        if arcpy.Describe("toselect").FIDSet: #If overlap = true, do the clip
            out_poly = os.path.join(outFolder, outName)
            arcpy.Clip_analysis("toselect", "ptmp", out_poly)
            arcpy.Delete_management("ptmp")
            arcpy.Delete_management("toselect")
            count = count + 1
        else:#Else, delete the tmp layers so they could be re-created for the next fc and the next row
            arcpy.Delete_management("ptmp")
            arcpy.Delete_management("toselect")
            count = count + 1 ## The script stop here


texts: TimeSlider on Web App does not match the TimeSilder on the underlying Web Map on ArcGIS Online

texts:
On ArcGIS Online, I have a hosted Web Map that contains 4 hosted feature layers. Each feature layer has the same number of features and a "Date" field and was published with Time enabled. All the features in a featuer layer has the same date and time for that field.

They have the following attribute tables.

The web map was published and configured to allow time slider to show the feature layers with time.

Then I created a web app using Web AppBuilder out of the web map with a TimeSlider widget added. The TimeSlider on the web app has three time steps while that for the web map has four steps as above. Missing the last time step stopped the last feature layer from appearing on the web app.

I thought the TimeSlider settings on a web app should match the settings on the underlying web map. Any idea?

Edit:
the URLs to them are:

The web map is http://arcg.is/1GHbaD

The web app made from hte web map is http://arcg.is/1vC5zm