texts: How to calculate great circle distance between unprojected SpatialPolygons?

texts:
I'm working (in R) with a SpatialPolygonsDataFrame describing US historical counties, and I would like to calculate great circle distances between centroids of counties (polygons).

I've looked into functions like spDists from the sp package, gDistance from rgeos, and distGeo from geosphere, but I'm confused about what exactly they're calculating. For instance, the documentation for spDists states that it calculates Euclidean distance if the spatial object is not projected, and great circle distance if the spatial object is projected. Similarly, gDistance rejects my data because it "expects planar coordinates". This is my first source of confusion because I would think it should be the opposite, eg isn't it more natural to calculate Euclidean distance from projected data because it is embedded in two-dimensional space? Similarly, shouldn't it be easier to calculate great circle distance from longitudes and latitudes (my data) than from coordinates for some projection that may have distorted distance?

I don't mind projecting my data to use these functions, but I'm not sure exactly what comprises a "projection." The output of proj4string for my data is

> proj4string(US)
[1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

I've seen in online tutorials that there's usually a "+init=..." component. Is that why it's not projected? Does it matter what "init" value I use to project it?


texts: Calculate real surface area with PostGIS

texts:
I would like to know if its possible to calculate a real surface area (with DEM information) only with PostGis function ?
I know it's possible with Saga and Grass.


texts: Python plugins don't work on QGIS 2.14 for Scientific Linux 7

texts:
I'm on Scientific Linux 7 and have been trying to install QGIS. According to the instructions, I need to

sudo wget fedora.vitu.ch/EL/qgis.repo -P /etc/yum.repos.d/
sudo rpm --import http://fedora.vitu.ch/Fedora/RPM-GPG-Key-vitu
sudo yum update
sudo yum install qgis qgis-python qgis-grass qgis-mapserver

The first step fails with a 403, but running the final two steps anyway did work, as QGIS seems to be available in the epel repo.

Now when I run QGIS though, I get the error "No module named sip".

I've looked everywhere and the only installers/source for SIP are for python 3. Is it possible to install SIP on python 2.7? Is there a workaround to get QGIS plugins working?


texts: Recommendations for road maintenance management (on a small scale)

texts:
I am in a road association and we maintain our own 3 mile dirt road.  In doing research on best practices for maintenance I discovered that other organizations use some mapping/GIS programs to help.

I know very little about GIS in general and wonder what would be good resources to start with.

More specifically, is there a generic program I can use to create a map of the road, then use it to collaborate with others to track issues with the road, map points on it and keep a history of maintenance?

If not, is there a base for this kind of thing and can a developer build one in a reasonable amount of time?  I am a software developer and would b more than happy to spend time on this project if that is required to get different pieces working.

The road is in Vermont.

EDIT:

Basically what I am hoping for is the following:

Software that allows me to drive a specific route and store that as a working model.
Download some GPS/County information to automagically set elevations, etc.
Allow multiple people to use the software and share information
Look at historical data over time (e.g. locations of potholes, bad ice, washboarding in a certain year or month - and view that as a report) - or some other means of collecting and storing data sets and keeping them distinct in time.

Please disabuse me of any unreasonable expectations or educate me on what I should be looking for.


texts:
I've got two datasets - 1) a set of buildings in LA and 2) fiber optic cable maps from a company called Zayo which covers the whole of the US.  I'd like to find the buildings within a certain distance to the fiber. The SRID is 4326 and it's stored as a geometry, not geography (as far as I can tell).

Here's my code:

select *
from la_records
inner join zayo_fiber
on ST_DWithin(zayo_fiber."geom", la_records."geom", .001);

The code looks to work, but  it takes forever - 10 minutes. Here's the Explain:

Nested Loop  (cost=0.21..43.81 rows=1 width=3167517)
->  Seq Scan on zayo_fiber  (cost=0.00..1.05 rows=5 width=3167299)
->  Index Scan using sidx_la_records_geom on la_records  (cost=0.21..8.54 rows=1 width=218)
    Index Cond: (geom && st_expand(zayo_fiber.geom, '0.001'::double precision))
    Filter: ((zayo_fiber.geom && st_expand(geom, '0.001'::double precision)) AND _st_dwithin(zayo_fiber.geom, geom, '0.001'::double precision))
Time: 0.045s

To speed up this query, should I break up the zayo fiber into smaller subsegments, then index it?


texts:
How to use FactoQGIS for Typological Analysis Principal Components Analysis (PCA) for vector features in QGIS?
I have a shapefile of polygons as in the following image:

Both visually and as illustrated by the Rose Diagram (in the image below) using Line Histogram Tool QGIS plugin available from the following link: https://plugins.qgis.org/plugins/LineDirectionHistogram/, it can be deduced that the overall orientation of the polygons is around 30 degrees, confirmed by calculating the mean of angle of orientation of all polygons.

I am looking for a Mathematical Tool to validate the mean angle.  I was advised to use the Principal Component Analysis (PCA) to find first eigenvector which gives the direction of the first principal direction.
I found FactoQGIS: a GUI tool based on an R script to perform Geometric Data Analysis in QGIS available from the following link: https://github.com/ESO-Rennes/FactoQGIS.
Below is the screenshot of FactoGIS Tool.

I have been trying the FactoQGIS Tool with different parameters, but could not generate results.
Which parameters should I set?
The shapefile of the polygons is available from the following link: https://drive.google.com/drive/folders/1BKp5VbR1aMAzCObGprvqeL0WSMqpTX32?usp=sharing


texts: Finding neighbouring countries using Overpass?

texts: How can I find all countries that are neighbours (share a border point or line) of a given country (eg. germany) using Overpass?


texts: Iso Cluster tool in script failing to open clip raster?

texts:
I have been coding in Python for only a few months, so maybe I'm just overlooking a minor detail as a newbie.

I have written the script below for a class project. It is supposed to clip two rasters, perform and unsupervised classification with the ISO tool, and then perform a change detection using map algebra. However, as you will see in error code below, the ISO tool is unable to open my clipped rasters. The script runs very fast. It's almost as if the ISO tool is attempting to open the rasters before the clip tool is done creating them. One two occasions, the script actually managed to produce the first classified image and only failed on creating the second one.

The script and error messages are included below:

## Import Necessary Modules

import arcpy
import os
from arcpy.sa import *

## Define repeat functions

def clip_it(AOI, raster, date):
     """This function clips file raster to supplied AOI file. It always uses
       AOI extent for clipping extent. Output file will be a TIF image with the
       supplied date as a filename"""
    try:
        arcpy.Clip_management(raster, '#', '{}clip.tif'.format(date), AOI, '0',     'ClippingGeometry')
    except:
        print('Clip Failed')
        print
        print arcpy.GetMessages()
        print
    return

def classify_it(raster, num_class, min_class, samp_int, date):
    try:
        outUnsupervised = IsoClusterUnsupervisedClassification(raster, num_class,     min_class, samp_int)
        outUnsupervised.save("{}class.tif".format(date))
        del outUnsupervised
    except:
        print('Classify failed')
        print
        print arcpy.GetMessages()
        print

    return

## Set up workspace

from arcpy import env
env.workspace = "u:/shared/gis/studata/jpdees0754/gitprojects/Final_Project/files"
env.overwriteOutput = True
path = env.workspace

## Prompt user for AOI and LANDSAT Rasters
## Promt user for raster data month and year
AOI = 'AOI.shp'
new_image = 'g2011Clip.img'
new_date = '2011'
old_image = 'g2000Clip.img'
old_date = '2000'
num_class = '5'
min_class_size = '20'
sample_int = '10'

## Clip both rasters to area of interest

clip_it(AOI, new_image, new_date)
clip_it(AOI, old_image, old_date)

## Perform unsupervised classification on raster images

classify_it('{}clip.tif'.format(new_date), num_class, min_class_size, sample_int,     new_date)
classify_it('{}clip.tif'.format(old_date), num_class, min_class_size, sample_int, old_date)

## Perfom change detection on classified rasters

new_class = arcpy.Raster("{}class.tif".format(new_date))
old_class = arcpy.Raster("{}class.tif".format(old_date))

outEqualTo = Raster(old_class) == Raster(new_class)
outEqualTo.save("{}_{}_Change.tif".format(old_date,new_date))

Error Message:

  Classify failed

  Executing: IsoClusterUnsupervisedClassification u:/Shared/GIS/StuData/jpdees0754/gitprojects/Final_Project/files\2000clip.tif 5 u:/Shared/GIS/StuData/jpdees0754/gitprojects/Final_Project/files\isocluster6 20 10 #
  Start Time: Tue Jul 29 22:10:23 2014
  Running script IsoClusterUnsupervisedClassification...
  ERROR 999999: Error executing function.
  Failed to open raster dataset
  Failed to execute (MLClassify).

  Error in executing Iso Cluster Unsupervised Classification.
  Completed script IsoClusterUnsupervisedClassification...
  Failed to execute (IsoClusterUnsupervisedClassification).
  Failed at Tue Jul 29 22:10:25 2014 (Elapsed Time: 2.00 seconds)


texts: Create closest point on multilinestring to use in shortest_path()

texts:
I am working a routing based application where I need to create a route between two points of interest. The problem which I am facing occurs after calculating network topology using assign_vertex_id().  When I run shortest_path I get the path between vertexes and not from the point of interest (as shown in image).

Can you please help me out with this: how can I use this query or any other function to do it?

Thanks for any tips.

EDIT 1

I am using this function called multiline_locate_point() similar to line_locate_point() which helps me to find nearest LINESTRING from MULTILINESTRING to the point.

here it is:

-- Function: multiline_locate_point(geometry, geometry)

-- DROP FUNCTION multiline_locate_point(geometry, geometry);

CREATE OR REPLACE FUNCTION multiline_locate_point(amultils geometry, apoint geometry)
  RETURNS geometry AS
$BODY$
DECLARE
    mindistance float8;
nearestlinestring geometry;
nearestpoint geometry;
i integer;

BEGIN
    mindistance := (distance(apoint,amultils)+100);
    FOR i IN 1 .. NumGeometries(amultils) LOOP
    if distance(apoint,GeometryN(amultils,i)) < mindistance THEN
        mindistance:=distance(apoint,GeometryN(amultils,i));
        nearestlinestring:=GeometryN(amultils,i);
    END IF;
END LOOP;

nearestpoint:=line_interpolate_point(nearestlinestring,line_locate_point(nearestlinestring,apoint));
    RETURN nearestpoint;
END;
$BODY$
  LANGUAGE plpgsql IMMUTABLE STRICT
  COST 100;
ALTER FUNCTION multiline_locate_point(geometry, geometry) OWNER TO postgres;

this gives me the_geom() of point on line. Now the problem is shortest_path() only takes vertex (integer number) as input to create route. So how can I pass this point to shortest_path() so that it creates route from the point and not from any vertex ??


texts: QGIS 1.8.0 Mac where is measure window or measure area symbol?

texts: Just trying to measure the area of a selected polygon (a few km2) but can't find out where the measure symbol is (mentioned in some posts here) or how to get the measure window shown in the manual. How?


texts: I am using ArcGIS 10.2.2 with the data interoperability extension. I am also building a web app using python, django and django-rest-framework. I'm building a json api that i can add to arcmap using data interop. but i can't figure out how to get arcmap/catalog to dump the cached json file it downloads and replace it with an updated version (as i change the api output). Any ideas?


texts: Copying one feature attribute to other feature using FME

texts: I have two features road and bridges. Road feature had an attribute UniqueID as string. I have a attribute field in bridge UponID. I want to copy road unique id to bridge Upon Id field. How to do it using FME?


texts: Multiple dimensions problem in importing kml files

texts:
When I import multiple kml files into R by using readOGR in a loop, I came across the problem denotes as follows

  Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
  use_iconv = use_iconv,  :  Multiple # dimensions: 2:3

Could anyone help me with this multiple dimensions problem?


texts: Is it possible to use GDAL to combine (rather than mosaic) rasters

texts:
Is there a GDAL command which can combine multiple rasters so that a unique output value is assigned to each unique combination of input values, like the combine command in ArcGIS?

gdal_merge seems to mosaic rasters so it is not exactly what I need.


texts: Remove features and dynamic vector layer from the map

texts:
Fixed the problem of showing or hiding the vector layer on the map after executing a query, I find that after a new search, the previous vector layer is not removed, it remains on the map.
I tried to make a vectorlayer.clear (), before adding the new features to the new vector layer, but when checking to show or hide the vector layer, at least the one generated when doing the previous search is always visible.
function drawAbundance(cbx_espId){
    var vectorAbundance = new ol.source.Vector({
        format: new  ol.format.GeoJSON(),
        loader: function(extent, resolution, projection) {
            var proj = projection.getCode();
            var url = "Abundance.php"
            var xhr = new XMLHttpRequest();
            var params = {
                abundance: cbx_espId,
            }
            xhr.open('POST', url);
            xhr.setRequestHeader("Content-Type", "application/json;charset=UTF-8");
            var onError = function() {
               vectorAbundance.removeLoadedExtent(extent);
            }
            xhr.onerror = onError;
            xhr.onload = function() {
                xhr.responseText;
                console.log(xhr.responseText);
                if (xhr.status == 200) {
                    **//I tried to put here:**
                    **//vectorAbundance.clear();**
                    vectorAbundance.addFeatures(
                        vectorAbundance.getFormat().readFeatures(xhr.responseText)
                    );
                    //map.getView().fit(vectorAbundance.getExtent());
                    $("#formEsp")[0].reset();
               } else {
                 onError();
               }
             }
             xhr.send(JSON.stringify(params));
        },
    });
    vectorlayer = new ol.layer.Vector({
        source:vectorAbundance,
        visible:true,
        displayInLayerSwitcher: false,
        style: style,
        maxZoom:9,
    });
    map.addLayer(vectorlayer)
}
//check to show or hidde the vector layer
function add_layer(evt){
    if(evt.checked){
        map.removeLayer(vectorlayer);
        //map.addLayer(vectorlayer);
    }else{
        //map.removeLayer(vectorlayer);
        map.addLayer(vectorlayer);
    }
}

Is there a way to make the previous search not keep on the map?


texts: shp2pgsql throws "ERROR: current transaction is aborted, commands ignored until end of transaction block"

texts:
I am trying to upload a shapefile to postgres using the command :

shp2pgsql -I -s 4326 "/Users/abc/Landmark.shp" table_name | psql -U user_name --dbname=db_name

I get the error :

  ERROR:  current transaction is aborted, commands ignored until end of
  transaction block
  ERROR:  current transaction is aborted, commands
  ignored until end of transaction block
  ROLLBACK
  ANALYZE

Postgis extension is added to the user postgres.
Not sure what is it that I am missing here.


texts: Clearing Python Console in QGIS using Python Command

texts:
I'm developing a python script in which I use several "print" statements to view the intermediate results in python console of QGIS. Each time after executing the program, I have to manually clear the python console using "Clear Console" option.

Is there a way out to clear the printed statements in python console using a command?

I tried using the following command from this answer

qgis.console.clearConsole()

But I am getting this error message

AttributeError: 'module' object has no attribute 'console'


texts: Classifying above ground LiDAR points in ArcMap?

texts:
I recently used the LAS Dataset toolbar to classify a .lasd file. The process classified the majority of points as unassigned with the rest being ground points. The study area has many tall buildings and tree, which did not receive the appropriate classification.

Is there a way of improving the quality of automatic classifications?


texts:
Setup: I have a project (UTM 84 zone 36N) with a digital elevation map (30 meters) of Israel and a GeoJSON loaded into QGIS. Both appear fine.

Task: I want to clip the DEM, using the GeoJSON. Because the GeoJSON is narrower than 30 meters, in places, I am using the warp tool, with the "CUTLINE_ALL_TOUCHED=TRUE" specification, in order to include all of the pixels.

Here is my code:

gdalwarp -wo CUTLINE_ALL_TOUCHED=TRUE -overwrite
-s_srs EPSG:32636 -t_srs EPSG:32636 -multi -dstnodata 0 -q -cutline
C:\...\geojsonSample.txt -dstalpha -of GTiff
C:\...\IsraelDEM.tif C:/.../clippedDEM.tif

Problem: sometimes, the resulting clip has appropriate elevation data, but does not appear on the canvas. Sometimes it gives me all 0s or infinite values for elevations. Any idea of what is going on?


texts: GPS and GIS for real time monitoring of an object

texts:
How is it possible to store coordinate of object in server in real time using feasible hardware ( can be used by student) GPS into database and retrived in real time??

What are the names of equipment ?
Any link to good resource?

The purpose of this experiment is to provide real time location of object through internet(web).


texts:
In QGIS 3.16, I have a layer of polygons corresponding to particular ecological habitat zones.  I'd like to perform a series of intersections in order to add information to each polygon as to where it is located, for example: what town it is in, what parcel, what catchment/drainage basin, what protected area (and several more).
The issue is that of course, many of my habitat polygons cross administrative and hydrological borders.  So if I intersect my layer1 with a layer containing town borders, the resulting intersection layer will duplicate each cross-border polygon, creating a separate entry for each town that the habitat polygon overlaps.
I'd like to avoid this duplication, or at least come up with a post-intersection trick to collapse the duplicates.  What I would want, ideally, would be a way to concatenate multiple towns (or basins, or protected areas) in the same field.  For example:

[Habitat name]-[Town]
Pond1         - Arlington, Lexington
Bog1          - Lincoln
Pond2         - Lincoln, Lexington

I've not found this trick yet.


texts: ArcPy-generated string formats badly in map text element

texts:
I've got a script that's pulling a bunch of attributes and concatenating the information into a long string.

text_description = "Address: {}\n\n".format(dictADDR[key])
text_description += "{} Customers.".format(len(listCustomer))
if countPole == 0:
    text_description += "No Poles\n"
else:
    text_description += "Poles:\n"
    text_description += "  Pole Quantity: {}\n".format(countPole)
(et cetera)

The goal is to include this string both in an emailed final report when the script finishes, and include it in a text element on a PDF map (displayed next to the parcel and an aerial map).

This works fine for the email, and gives something like:

Address:  1225 BICKLEY RD
0 Customers
No Poles
No Overhead Distribution

However, in the text element it ends up ignoring spacing and overlaps some text with other text (see image below).

Copy-pasting the contents of the text properties box, it appears to have the line breaks in the correct places:

Address:  1225 BICKLEY RD
0 Customers.
No Poles
No Electric Overhead Distribution

What's going wrong, and what should I change in either the text element or my code to fix it?


texts: Converting map coordinates using AutoCad to Excel sheet

texts:
I have some maps of pieces of lands and I want to convert the coordinates to excel sheet

I want to arrange the data like this in the excel sheet:

land number( the pink number in the pic)  accompanied by its coordinates  in each row

I will be really thankful if you can help me

link to the map file above:

update

I want to have 3 columns in the excel sheet

1) representing the land number text value

2) representing the the coordinates of that land separated by "_" like this
 23.87718162405882_57.10554858573619_23.877357092058385_57.105595188322354_23.877292855536815_57.10588146117277 each two representing a point/corner

3) the block number where the land is

if the format of the coordinates can not be like above then there is no harm having 3 columns x and y or easting and northing plus point number. therefore each land will have rows of the same points number


texts: Seeking VIIRS Day Night Band and alternatives for night satellite images

texts:
Has anyone worked with VIIRS Day Night Band?

I want to download few
images to get a map of night scene of Central Europe (city lights etc).
I use friendly interface at NOAA (http://www.class.ncdc.noaa.gov/saa/products/search?datatype_family=VIIRS) and I've downloaded few images already. The problem is that there was only one nice image and the rest of them were in 95% covered with
value 0 (although they were taken in the same night!).

My questions are:

Is there any way to preview VIIRS Day Night Band images before download?
Does mosaic of night scene images for Europe already exist and is available for download?
Do you have any other ideas where to find remotely sensed images of earth at night?


texts: Calculating distance of points along route using ArcGIS Desktop?

texts:
I am struggling with the calculation of the distance between points along a stream (points are located 0-20m away from the stream) in ArcGIS 10. The points are simple xy-data in UTM.

I created a route from the shapefile of the stream and used the tool "locate features along routes" to add the points but how can I now calculate the distances of the points along the route?

Or is there an other tool how to do such a calculation?


texts: Finding polygons of one layer that either fit inside or overlap another specific polygon using PyQGIS

texts:
I have a layer with the population density of a city, in polygons that represent the city blocks. I also have a layer of arbitrary hexagons. I want to calculate the population that each hexagon has by finding the polygons that overlap each hexagon, iterating on those to find the area that is actually inside the hexagon and multiplying that by the population density.
First, I am trying to count the blocks that intersect a hexagon, just to get this working. But the QgsSpatialIndex.intersects(geometry) function returns a null result. This is my code:
from qgis.core import QgsProject
from qgis.utils import iface
from PyQt5.QtCore import QVariant

_Layer_hexa = 'Hexagon_Layer'
_Layer_densidad = 'Population_density'

layer_dens = QgsProject.instance().mapLayersByName(_Layer_densidad)[0]
layer_hexa = QgsProject.instance().mapLayersByName(_Layer_hexa)[0]

_DENS_FIELD = 'Dens pob'
_NEW_SUM_FIELD1 = 'Blocks'

# Create new field:
layer_hexa.startEditing()
layer_hexa.dataProvider().addAttributes(
        [QgsField(_NEW_SUM_FIELD1, QVariant.Int),
         ])
layer_hexa.updateFields()

# Feature dictionaries
feature_dict_hexa = {f.id(): f for f in layer_hexa.getFeatures()}
feature_dict_dens = {f.id(): f for f in layer_dens.getFeatures()}

# Build spatial indexes
index_hexa = QgsSpatialIndex()
for f in feature_dict_hexa.values():
    index_hexa.insertFeature(f)

index_dens = QgsSpatialIndex()
for f in feature_dict_dens.values():
    index_dens.insertFeature(f)

# Iterate the hexagons:
c=0
mzas = 0
for f in feature_dict_hexa.values():

    # Get the hexagon geometry
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.

# This is where it returns null------------------------------------------------------------
    intersecting_ids_dens = index_dens.intersects(geom.boundingBox())
#------------------------------------------------------------------------------------------
    #Initialize the sum:
    mzas = 0
    # Iterate on the intersections:
    for intersecting_id in intersecting_ids_dens:
        #Find the feature in the dictionary
        intersecting_f = feature_dict_dens[intersecting_id]

        # Check that the geometries actually intersect:
        if not intersecting_f.geometry().disjoint(geom):
            mzas += 1

    f[_NEW_SUM_FIELD1] = mzas
    c += 1
    # Updateamos la layer con este cambio
    layer_hexa.updateFeature(f)

layer_hexa.commitChanges()
print('Processing complete.')

I have read the QgsSpatialIndex documentation but am still at a loss as to why it isn't working. While this gets posted I will try to merge the hexagon and density layer and see if this works using only one layer, but I would still like to understand why the code above doesn't work.
I am basing my code around this code by Ujaval Gandhi for finding neighboring polygons: https://www.qgistutorials.com/en/docs/find_neighbor_polygons.html


texts:
I'm using Mapserver 7.4.0-dev and MapScript PHP.

In previous versions of MapServer and MapScript (https://mapserver.org/mapscript/php/index-5.6.html#shapeobj), I could query a shapeObj after its values. And get a key => value array with all the data fields.

Example:

public function __construct(\shapeObj $shape)
{
    foreach ($shape->values as $key => $value) {
        //Do something;
    }
}

In the new Mapserver and MapScript (https://mapserver.org/mapscript/mapscript.html#shapeobj), there is no $shape->values attribute.

I can get the value of element X but no the key name (Example "Gis_ID") by using this function:

$shape->getValue(3);

But this does not help my much. My map sources are dynamic and i need to target the data by IDs, and not a random array index.

Can anyone help me figure this out? :)


texts: Convert points to lines with Python GDAL


texts: What is the best way to handle geometries between Python functions?

texts:
What is the best way to handle geometries between different Python functions?
For example:
I have a function that exports a line from grass.script. The line will be a parameter of another function where it will be passed in.

Is it best to export them to WKT, or safe them in a GPX, or GeoJSON? The formats OGR offers are endless. But which one is most recommended?


texts:
I've installed pyModis on my Mac using pip install pyModis. No error messages. When I try to run modis_download.py I get this error:

Traceback (most recent call last):
  File "/usr/local/bin/modis_download.py", line 30, in <module>
    from pymodis import optparse_required
  File "/Library/Python/2.7/site-packages/pymodis/__init__.py", line 3, in <module>
    from . import downmodis
  File "/Library/Python/2.7/site-packages/pymodis/downmodis.py", line 53, in <module>
    from urllib.request import urlopen
ImportError: No module named request

Does anyone know what's up here?


texts: How to compare two rasters and reclass the greater value(python script)?

texts:
I am working on an assignment that requires us to reclass a raster(via python) from two rasters. We need to compare the power received at one cell(on raster1) and compare it to the same cell on raster 2.

If the value of that cell in raster1 is greater than the value in raster2, the cell would be assigned a value of "1" while a value in raster2 greater than raster1 would get a "2".  We also need to reclassify certain values as "0" since they are considered out of range.

So the final raster would have three values:
0-out of range/no service
1-better value in raster1
2-better value in raster2

I can try to clarify further if needed.


texts: Creating centerline for rasterized polygons using ArcScan?

texts:
I would like to use the automatic vectorization feature in ArcScan to create a centerline (line feature) for each of 20 rasterized polygon features.

Using a bi-color raster (created with polygon to raster tool and 1= foreground, 0= background values), and with centerline vectorization method selected from the ArcScan toolbar, I generated features from the vectorization tab on the ArcScan toolbar (and stored the output in a blank line feature class). However, the output created 3 small line features for only the smallest regions of the rasterized polygons and failed to create centerlines running the full length of my rasterized polygons.

I have read the related Creating Centrelines from Road Polygons/Casings? but did not find sufficient detail to create centerlines with ArcScan

I am using ArcMap 10.4.1 for desktop with an Advanced license and Spatial Analyst and ArcScan Extensions.


texts: Carto.JS - Pop up closes when mouse hovers out of feature

texts:
As per the CARTO docs I am using the following code to trigger a info pop up when the user's mouse hovers over a feature:

const popup = L.popup({ closeButton: false });
    ecosystem.on(carto.layer.events.FEATURE_OVER, featureEvent => {
      popup.setLatLng(featureEvent.latLng);
      if (!popup.isOpen()) {
        popup.setContent(featureEvent.data.name_2);
        popup.openOn(map);
      }
    });

However, the info pop up doesn't close when the user's mouse is moved out of the feature's bounds. The user has to click for the windows to close. Is it possible for the info pop up to only display when the mouse hovers over the feature, and closes automatically?

Here is a working example.


texts: Error [37] : Unable to calculate. Metadata error - QGIS SCP plugin

texts:
I was trying to preprocess the data I downloaded from Sentinel 2 through SCP plugin but Error 37 Shows up. I don't know why. I was suppose to clip the jp2 images but it also doesn't work readily so I guess I have to preprocess it first to get a certain type of file that the clip through SCP plugin can work with...
I'm Using Macbok - EL capitan
QGIS 3.8


texts: How to make autoIncrement number start with (001....009, 011......099)?

texts:
i have a challenge with my time in my work, i have feature class called "Parcel" storing many data, i generated an autoIncrement number by using field calculator, in field "parcel code"

Expression:
autoIncrement()

Expression Type:
PYTHON_9.3

Code Block:
rec=0
def autoIncrement():
 global rec
 pStart = 1 #adjust start value, if req'd
 pInterval = 1 #adjust interval value, if req'd
 if (rec == 0):
  rec = pStart
 else:
  rec = rec + pInterval
 return rec

that i asked here before
Adding automated values to attribute field?

i want to start number to 001......009, and 010.........099,


texts:
I have a simple table (location_test) with three column (name, geom , jid). geom is a geometry column with SRID 27700. I set the user's geometry position with:

UPDATE location_test SET geom = ST_SetSRID(st_makepoint(48.9 , 2.4),27700) WHERE location_test.jid = 'luca';
UPDATE location_test SET geom = ST_SetSRID(st_makepoint(48.9 , 2.4),27700) WHERE location_test.jid = 'filo' ;
UPDATE location_test SET geom = ST_SetSRID(st_makepoint(48.9 , 2.4),27700) WHERE location_test.jid = 'gianlu';

UPDATE location_test SET geom = ST_SetSRID(st_makepoint(55.8 , 37.6),27700) WHERE location_test.jid = 'mango';

3 people are in Paris (lat 48.9 ,lon 2.4) and the last one is in Moscow (lat 55.8 ,lon 37.6).

When I query for people in Paris, within 100 meters, "mango" is included, but why? the query is:

SELECT * FROM location_test WHERE ST_DWithin(geom, ST_GeomFromText('POINT( 48.9  2.4 )',27700) , 100.0);

in SRID 27700 units are meters, so 100.0 is in meters i think.

could someone help me?


texts: How to make query task synchronous in ArcGIS JSAPI

texts:
I am using ArcGIS JS API version 3.6. I want to make query task synchronize. While executing query task reaming code will wait until completion of query task.

For example I have shared my code here which contains query layer.

Case 1: In which I am directly executing the code and printing the log in console window, where you can see before the query task output is displaying

Case 2: In this case we have used alert (for time delay) in the code then value is getting before query task.

-Researching  on Google and going through similar question

Case 1 screenshot :

Case 2 screenshot :

Any help will be great !!! thanks in advance : )


texts:
So I happen to have a little bit of web experience, and have recently bought an OBD-II device. Basically, I'm trying to log certain metrics from the car, such as location and other diagnostics.

I can send in an SMS to the SIM on the device to configure the server IP that the device sends data to... But I'm not sure what happens next.

A piece of java code that interprets info from the device and makes it available to the PHP & JavaScript applications?

I'm using this: eeLink OBD - GPS Device...


texts: Create a way to easily select by attributes for a complicated union file

texts: I have polygon feature class created from a union of 3 other feature classes.   These are 5 ownership types, 3 classes of slope % and a multi-ring (3 distances; 50,75 and 100) road buffer.   So any given polygon could have 1 of 45 combinations of features; 3*5*3 =45.    I'd like to make a simple tool where a user could pick the combinations of features to quickly kick out only those areas as a mask to be used in another program.


texts: Web Cartography how to, CartoDB - Mapbox

texts: I have to prepare an online map like the following http://www.oikoskopio.gr/map/ but I don't know what program they use. Can I do something like that with CartoDB or Mapbox; with these services you can make a map but the TOC, can I do something like this in CartoDB or in Mapbox?


texts: How to split Indian states and provinces from world states and provinces shapefile?

texts:
I am following Mike Bostock's D3 Lets make a map to create GeoJSON for states and provinces of a country. I got a shapefile for states and provinces from Natural Earth. I am trying to extract info for India with this command:

ogr2ogr -f GeoJSON -where "ADM1_A3 IN ('IND')" subunits.json ne_10m_admin_1_states_provinces.shp

Here is the error:

  GeoJSON driver failed to create subunits.json

Is there any ISO country code for Admin-1 level?


texts: Editing Labels in Legend using arcpy in pycharm

texts:
I have created this map in PyCharm with Arcpy software, and was wondering if there is a code I can use to remove the 'High' and 'Low' labels and just show the numeric values in the legend. Using ArcMap version 10.7.1.
Here is a sample of the code I've used to edit items in the legend:
#Using Style Manager to edit the appearance of the Legend
#edit study area in legend
mxd1 = arcpy.mapping.MapDocument(r"C:\\Users\\Windows\\Desktop\\test\\map\\map_test.mxd")
df1 = arcpy.mapping.ListDataFrames(mxd1, "Layers")[0]
lyr1 = arcpy.mapping.ListLayers(mxd1, "Noise (dBA)", df1)[0]
styleItem = arcpy.mapping.ListStyleItems("C:\\Users\\Windows\\AppData\\Roaming\\ESRI\\Desktop10.7\\ArcMap\\Windows.style", "Legend Items", "Noise (dBA)")[0]
legend = arcpy.mapping.ListLayoutElements(mxd1, "LEGEND_ELEMENT", "Legend")[0]
legend.updateItem(lyr1, styleItem)
mxd1.save()

I have also used the 'apply symbology from layer' tool, but this did not effect or change the labels.


texts: Determining proportion of raster cells containing each value from list for many rasters in raster stack?

texts:
I wrote a function to determine the proportion of raster pixels that contain each of a list of values that I would be interested in.  I then wrote a function to apply the previous function to a 3d numpy array.  I would like to know how to improve this and also whether it would be easier to apply a function like this to a list of raster files, rather than create a numpy array.

A test array and sample values:

import numpy as np
test_arr = np.random.randint(0, 200, 200).reshape(2,10,10)
values = [test_arr[1][1][1], test_arr[0][0][0]]

Here are my functions:

def pixel_props(rast, values):
"""this function outputs a list with the proportions of
the total raster pixels that contain each value."""
prop_list = []
size = float(rast.size)
for i in values:
    temp_count = np.count_nonzero(rast[rast == i])
    prop_list.append(temp_count/size)
return prop_list

def raster_props(array, no_data):
"""this will apply the pixel_props function to a bunch of layers in a numpy 3d array"""
    master_matrix = [[]]
    master_matrix.append(no_data)
    for layer in array:
        master_matrix.append(pixel_props(layer, no_data))
    return master_matrix


texts: How to transform y-axis values of Google Earth Engine chart

texts:
I'm using Chart.image.series() in Google Earth Engine to create a line chart with two series from an image collection of processed land cover data (two binary bands per year: grass/shrub/wetland, and cropland/developed) and a region (10 mi radius circle).

I would like the y-axis values (sum area) to be log transformed but I can't figure out how to do that. Such as if ee.Reducer.sum().log() was allowed.

https://code.earthengine.google.com/c08b3d0316375648af0fe1cd4fe8128a

//try to calculate landcover rate of change for cropland/developed and grass/shrub/wetlands

//Get CDL, and point of interest
var imageCollection = ee.ImageCollection("USDA/NASS/CDL")
var geometry = /* color: #d63000 */ee.Geometry.Point([-100.12247612030984, 46.13719639996092]);
var mi10 = geometry.buffer(16093.4);

Map.centerObject(mi10, 10);

var ic= imageCollection.filterDate('2008', '2019');
//print(ic);

//get grass, shrub, wetlands
//other hay non-alfalfa 37, shrubland 64, 152, wetlands 87,  190, 195, grass 176
//get crop and development
// crop 1-36, 38-61, 66-77, 82, 204-254, devleoped 121-124
var gogetit = function(image){
  var cropland = image.select('cropland');

  var total = cropland.eq(37).or(cropland.eq(64)).or(cropland.eq(152))
  .or(cropland.eq(87)).or(cropland.eq(195)).or(cropland.eq(190)).or(cropland.eq(176));

  var cropdev = cropland.gte(1).and(cropland.lte(36))
  .or(cropland.gte(38).and(cropland.lte(61)))
  .or(cropland.gte(66).and(cropland.lte(77))).or(cropland.eq(82))
  .or(cropland.gte(121).and(cropland.lte(124)))
  .or(cropland.gte(204).and(cropland.lte(254)));

  var total_fm = total.focal_mode(75, 'square', 'meters');
  var cropdev_fm = cropdev.focal_mode(75, 'square', 'meters');

  return total_fm.addBands(cropdev_fm).rename(['grass/shrub/wetland', 'crop/developed'])
  .set('system:time_start', image.get('system:time_start'));
};
var ic2 = ic.map(gogetit);
print(ic2)

var crop08 = ic2.first().select('crop/developed');
var grass18 = ee.Image(ic2.filterMetadata('system:index', 'equals', '2018').first());
var grass18 = grass18.select('grass/shrub/wetland');
var grass08 = ic2.first().select('grass/shrub/wetland');
Map.addLayer(crop08.updateMask(crop08),{palette:['yellow']}, 'crop08')
Map.addLayer(grass08.updateMask(grass08),{palette:['red']}, 'grass08')
Map.addLayer(grass18.updateMask(grass18),{palette:['green']}, 'grass18')
Map.addLayer(mi10, {}, 'area');

//convert to area
var prep = function(image){
return image.multiply(ee.Image.pixelArea()).set('system:time_start', image.get('system:time_start'));
};
var ic2 = ic2.map(prep);
//print(ic2);

//line chart for two landcovers with trendlines
var chartit = new ui.Chart.image.series(ic2, mi10, ee.Reducer.sum(), 30)
.setOptions({
      title: 'Grass/Shrub/Wetland vs. Crop/Developed',
      hAxis: {'title': 'Year'},
      vAxis: {'title': 'Area'},
      trendlines: { 0: {showR2: true, visibleInLegend: true} ,
                    1: {showR2: true, visibleInLegend: true}}
                    });
print(chartit);


texts: Validate coordinates in DMS format in a Textbox

texts: I have to validate a textbox which is taking coordinates in Degree Minutes and Seconds or Degree and Minutes or in Degrees format divided by a space between them.


texts: Raster functions performing poorly in ArcGIS Desktop?

texts:
After installing 10.5 I noticed some documents were taking much longer to open, some 4 to 15 minutes and others would not open at all.
Initially I thought some databases/rasters were corrupted but it took some troubleshooting and a machine that still runs 10.4 to realize that the problem is in 10.5 specifically any document that contains a raster function layer even if it is turned off takes exponentially longer to open / refresh.  (the more raster function layers in the document, the worse the problem)

Removing all the raster function layers solves the issue and the documents suddenly become as snappy as before.  Naturally, the next step was to attempt to rebuild these raster functions from scratch in 10.5.  Using the Image Analysis window and simply selecting a hillshade and an imagery raster and hitting the pansharpen buton to create a raster function layer takes MUCH longer in 10.5. It does work and eventually a raster function layer is loaded in the TOC, however, is take a very long time and further attempts to tweak the function / raster function chain is not possible because the Function tab in the Layer Properties crashes.

Has anyone encounter similar issues and if so, is there a workaround? Because most of my 100s of documents contain raster functions, I may have to downgrade to 10.4.1 but I am hoping for a more elegant solution. I already created and resaved a number of documents which makes these incompatible with legacy ArcMap versions so I will need to resave all these as 10.4 copies first.

Update: After a back andf forth with ESR, this was confirmed and is now a BUG-000103911.  The workaround is: Downgrade to 10.4.1


texts: PostGIS - Linestring doesn't contain a Point - Problem

texts:
I have a problem to check if a Linestring contains a Point.
I consider the straight line y = 0.1 x , and take by it a segment l='Linestring( 0 0 ,  4 0.4 )' and a Point p on segment l, p = 'Point(1.5 0.15)'. If I check ST_Contains(l,p) the result is false and the distance  ST_Distance(l,p) is near to zero, but this is wrong because the Point is on the Linestring.
If I try with p= 'Point(2 0.2)' the result is right.

What is the problem?

WITH
     l AS (
           SELECT ST_GeomFromText('Linestring( 0 0 ,  4 0.4 )') AS geom          -- segment of the straight line y = 0.1 x
       ),
     p AS (
       SELECT ST_GeomFromText('Point(1.5 0.15)') AS geom                     -- point on the straight line y = 0.1 x
       )

 SELECT ST_AsText(ST_Intersection(p.geom,l.geom)) as ST_Intersection,
        ST_Distance(p.geom,l.geom),
        ST_Contains(l.geom,p.geom),
        ST_Intersects(p.geom,l.geom),
        ST_AsText(ST_Line_Interpolate_Point(l.geom,1.5/4)) as ST_Substring       -- point in position 1.5/4 of l
 FROM l,p

EDIT:
My origin problem was:
LineString L1 contains Linestring L2 but if I do ST_SymDifference(l1.geom,l2.geom) the result is incorrect.

WITH
     l1 AS (
           SELECT ST_GeomFromText('Linestring( 0 0 ,  4 0.4 )') AS geom          -- segment of the straight line y = 0.1 x
       ),
     l2 AS (
       SELECT ST_GeomFromText('LineString(0 0 , 1.5 0.15)') AS geom             -- segment of the straight line y = 0.1 x
       )

 SELECT ST_AsText(ST_SymDifference(l1.geom,l2.geom)) as ST_Intersection
 FROM l1,l2


texts: Passing File Directly to esriRequest Instead of Using Form Upload?

texts:
The esriRequest is used for retrieving data from a remote server or uploading a file. The request has a Form option which is used if the request is to upload a file. We specify the form element that contains the file input control here (as shown below).

request({
    url: portalUrl + '/sharing/rest/content/features/generate',
    content: myContent,
    form: dom.byId('uploadForm'),
    handleAs: 'json',
    load: lang.hitch(this, function (response) {
    if (response.error) {
      errorHandler(response.error);
      return;
    }
    var layerName = response.featureCollection.layers[0].layerDefinition.name;
      addShapefileToMap(response.featureCollection);
    }),
    error: lang.hitch(this, errorHandler)
});

Now, I need to pass a file directly to the options instead of using HTML File Dialog and Input. Something like:

var data = "www.domain.com/GIS/App.ZIP";
 request({
 ....,
 form: data,
 ....
});


texts: What preprocessing is required when using Landsat for LULC?

texts: Can anyone help me regarding the what kind of image pre-processing techniques should be used for the landsat images downloaded from http://earthexplorer.usgs.gov/ to create land-use land-cover (LULC) map?


texts: Extracting only part of supervised output model in Raster format

texts:
I have a svm predicted model which when applied on raster images classifies the LULC of land and writes the output to raster file. I want to use only the Land used part out of the output raster image. Any leads on this.
Below is the code used for svm prediction:
bands_data = np.dstack(bands_data)
bands_data = np.nan_to_num(bands_data)
rows, cols, n_bands = bands_data.shape
n_samples = rows * cols
flat_pixels = bands_data.reshape((n_samples, n_bands))
r = clf.predict(flat_pixels)
result = np.array([int(i) for i in r])
classification = result.reshape((rows,cols))

def write_geotiff(fname, data, geo_transform, projection):
    """Create a GeoTIFF file with the given data."""
    driver = gdal.GetDriverByName('GTiff')
    rows, cols = data.shape
    dataset = driver.Create(fname, cols, rows, 1, gdal.GDT_Float32)
    dataset.SetProjection(projection)
    band = dataset.GetRasterBand(1)
    band.WriteArray(data[np.where(data == 1)])
    dataset = None #Close the file
    print("done")

output_file = fname
write_geotiff(output_file,classification, gt, proj)


texts: Custom osm2po postprocessor to generate a node file and a way file

texts:
For a software component were using we need to provide 2 input files. One with a list of nodes (nodeId, lat, lon) and one with ways (fromNodeId, toNodeId, oneway, maxSpeed, lenght).

To create these 2 files I created a custom postprocessor class, you can find the code here: https://gist.github.com/anonymous/5047904

To get the Nodes I read the "tn_raw.2po" file, for the ways the "sw_all.2po" file. This works great, but now my node output files contains all osm nodes. What I would like, is that the file only contains the nodes used in the ways file.

Is it possible to only output the used nodes in a file?

Thanks!


texts: CartoDB st_makeline with geography datatype

texts:
I am using an online cartodb account to try and create a shapefile based on flight routes between one airport and another.

So far, I have uploaded my Lat/Long points for each connected airport, and used st_makeline to create a line string between them in geometry data type.

However, this outputs a direct straight line, whereas it should really follow a geographic projection and output as a curved line.

Does anyone have any advice on how to do this? Bearing in mind I am using cartodb online so am restricted to the custom SQL query it provides.

Current query is as follows:

SELECT
    routeid,
    boardpoint_offpoint,
    ST_Makeline(the_geom) as the_geom
FROM
    routemapping
group by
    routeid, boardpoint_offpoint


texts: Increase rendering speed for the data coming from MSSQL for showing in Leaflet

texts:
I have some polygon data which are coming from MSSQL server. Getting these data takes about 5 seconds and I have about 1600 Polygons. I would like to show a popup for each polygon (by clicking on them, helper method). The rendering time after getting this data for showing them in a Leaflet map takes about 40 seconds which is too long. How I can accelerate the rendering process?

My code:

What I have tried:

Change basemap providers, I used OSM, not a dramatic improvment.


texts:
I added the popup (click) on the data layer, now how can I add the hover popup without broke the click popup?

The hover popup it will be just a few thing, not the whole text I display with the click popup, it will be like a preview.

Here my JS code:

    function info() {

      map.on('click', function(evt) {
        var content = document.getElementById('popup-content');
        content.innerHTML = '';
        var viewResolution = map.getView().getResolution();
        var coordinate = evt.coordinate;
        url = pmfeatlayer.getSource().getGetFeatureInfoUrl(evt.coordinate, viewResolution, projection, {'INFO_FORMAT': 'text/html'});

        if (url) {
          var xhttp = new XMLHttpRequest();
          xhttp.open("GET", url, true);
          xhttp.send();
          xhttp.onreadystatechange = function (aEvt) {
            if (xhttp.readyState == 4) {
              if (xhttp.status == 200)  {

                var parser = new window.DOMParser();
                var res = parser.parseFromString(xhttp.responseText, "text/xml");

                if (window.innerWidth < 768) {
                  // show hide modal on mobile and if consult or not
                  if(res.getElementsByTagName('item').length != 0){
                    if (window.location.href.indexOf("consult") != -1) {
                      fullInfoPc(xhttp.responseText);
                    } else {
                      fullInfo(xhttp.responseText);
                    }
                  } else {
                    $("#mobpopup").modal('hide');
                  }

                } else if (window.innerWidth >= 768) {
                  // set the position of the popup on the screen (from tablet) and if consult or not
                  if(res.getElementsByTagName('item').length != 0){
                    overlay.setPosition(coordinate);
                  } else {
                    overlay.setPosition(undefined);
                  }

                  // create the content of the popup
                  if (window.location.href.indexOf("consult") != -1) {
                    fullInfoPc(xhttp.responseText);
                  } else {
                    fullInfo(xhttp.responseText);
                  }

                } // end else if
              } // end if 200
            } // end if 4
          }; // end onreadystatechange

        } // end if(url)
      }); // end map.on(click)

      map.on('pointermove', function(evt) {
        if (evt.dragging) {
          return;
        }
        var pixel = map.getEventPixel(evt.originalEvent);
        var hit = map.forEachLayerAtPixel(pixel, function() { return true; }, null, function(layer) { return layer === pmfeatlayer; });
        map.getTargetElement().style.cursor = hit ? 'pointer' : '';
      });

I'm trying to use this example but I can't make it works, someone know how to do it?


texts: Access feature data for all features in a NamedMap

texts:
For any Cartodb NamedMap, we load it using this function:

  cartodb.createLayer(map, namedLayerSource, {https: true})
        .addTo(map)
        .done(function (layer) {

});

After that I need to access all features data inside that named map.
So It is known that we can do this:

sublayer.on('featureClick', function (e, pos, latlng, data) {
                  console.log(data);
                    }
                );

But what can i do to get data for all features without "featureClick" event ?
Is there anyway I can get all data of that NamedMap in an Array or a Map?

I am using Cartodb.js v3.15 And I want to enable a search feature in which the user can search for any details inside the data for all points/features exists in a carto NameMap.


texts: OpenLayers - draw feature with style, lost after map redraw

texts:
The function that I am trying to work on is: click on a feature and change its color. I tried to approach this by first recreating the feature with custom style:

layer.drawFeature(feature, myStyle);

this gave me the desired shape but it was lost after layer redraw.

Then I tried to use a clone of the feature:

layer.drawFeature(feature.clone(), myStyle);

this time the new shape would not disappear but it's not sync with the map, i.e. it does not resize/ move with the map.

I do notice that the following post was about a similar issue. I wonder if there's a way without creating an extra overlay? I could live with the cloned feature thing if that can be put to work.

Openlayers - Losing selection or style on layer redraw

Any suggestions will be greatly appreciated!


texts: QGIS, Calculate networks first, diameters second

texts:
I need to get the following task done:
First, various points (stored in a shapefile, identified by a common attribute) need to be put into relation with each other (see graphical scheme, steps 1 and 2). Ideally, this would result in a circle encompassing all points. This is where I am a bit lost at the moment because Thiessen/Voronoi/Delaunay polygons aren't of much help here.
Second, and this seems rather easy, the diameter of the resulting polygon should be written out in a text file (see steps 3 and 4). The main goal, as you may imagine by now, is to come up with some idea on what the average diameter of a network is. Therefore, if there is any better way to answer this, please don't hesitate to recommend it to me.
I am not familiar with too much Python. However, since this kind of analysis would need to be done a couple of hundred times, I am more than willing to spend some time on learning some scripting/automation.


texts:
What’s their relationship or difference?
I’ve read that WGS84 uses ECEF, but they can be converted between one another. If not mistaken, they are both geographic coordinate systems and not projected coordinate systems. But that’s about all I understand.


texts:
I was generating several hillshade effects from a small Projected Raster (WGS 84 - Zone 23S)and one of those failed. The script that failed was the Sky Model from the Terrain Tools 1.1 package (https://www.esri.com/arcgis-blog/products/3d-gis/3d-gis/terrain-tools-1-1-released/).This specific script runs based in a a pre-fabricated command line avaible with the package.
I ran a Sky Model using the Overcast script that came with the toolbox and receive this error message:
Traceback (most recent call last):
  File "D:\USUARIOS\Andre\Downloads\TerrainTools_v1_1\Scripts\SkyModel.py", line 81, in <module>
    skymodel(*args)
  File "D:\USUARIOS\Andre\Downloads\TerrainTools_v1_1\Scripts\SkyModel.py", line 63, in skymodel
    prevRaster.save(outRaster)
RuntimeError: ERROR 999998: Unexpected Error.
Failed to execute (SkyModel)

So, I tried several solutions and utimately format and reinstalled the Windows 10 and the ArcGIS 10.5. Didn't worked!
How can I fix it?


texts: Set CartoCSS ranges dynamically for choropleth

texts:
I know that I can generate a choropleth using CartoCSS like this:

#volume [ number_of_cases <= 434] {
   marker-fill: #BD0026;
}
#volume [ number_of_cases <= 183] {
   marker-fill: #F03B20;
}
#volume [ number_of_cases <= 130] {
   marker-fill: #FD8D3C;
}
#volume [ number_of_cases <= 90] {
   marker-fill: #FECC5C;
}
#volume [ number_of_cases <= 50] {
   marker-fill: #FFFFB2;
}

That works fine. However, I'm changing the field that sets the volume dynamically. For example, it might be the number of "widgets" and then it might be the number of "thingies". What I can't have is these ranges in the CartoCSS set the same for each product (widgets and thingies). The ranges need to be set "on the fly" based on the possible ranges in the number_of_cases column, as they have very different ranges. What I need is something that generates the ranges similar to how the CartoDB editor works, but to do it all dynamically.


texts: How can I count the number of countries my GPS data fall into?

texts:
I'm working with GPS data in R from different studies and would like to list the countries that each study covers. I can easily map the data using leaflet but I'm not sure how to extract the countries.

Hope you can help

#' load in the data

mydata <- mydata %>%
  convert(num(long:lat))
str(mydata)

#' map it

 m <- leaflet(mydata) %>%
  addTiles()  %>%
  # setView( lat=-27, lng=170 , zoom=4) %>%
  addProviderTiles("Esri.WorldImagery") %>%
  addCircleMarkers(
    ~ long,
    ~ lat,
    fillOpacity = 0.7,
    color = ~ pal(species),
    radius = 3,
    stroke = FALSE
  ) %>%
  addLegend(
    "bottomright",
    pal = pal,
    values = ~ species,
    title = "Species",
    opacity = 1
  )
m


texts: Performing ellipsoid to spheroid conversion using ArcMap

texts:
I have a DTM that is referenced to an ellipsoid (equatorial radius is not equal to the polar radius) and I want to convert the elevation values to the spheroid (equatorial radius = polar radius).
How can this be done using ArcMap?
Additional info: The DTM is of a different planetary body.


texts: ArcPy script doesn't work for TRIM function

texts:
I have this code:

import arcpy
import os

def massiveTrim (vrb):
    y=r"shape path"
    for i in y:
        campos = arcpy.ListFields(vrb, "", "String" )
        for camp in campos:
            espacio = "Trim( ["+str(camp.name)+"] )"
            arcpy.CalculateField_management(vrb, str(camp.name), espacio, "VB", "")

print "Finish process"

The script use the TRIM function for delete before and after spaces in a text string for all fields in a shape or FC, this task is automatic, the problem is that TRIM doesn't work.
I have ArcGIS 10.5, Pyscripter


texts: Merging multiple shapefiles into one NULLs out values

texts:
I am working in QGIS 3.4, and trying to merge 10 polygon shapefiles into a single shapefile.  However, the result is not preserving the values for most of the shapefiles, it's NULL-ing them out for all but one of the shapefiles.

Must the exact number of fields exist for all 10 shapefiles?
If so, must the fields be ordered exactly the same for all 10 shapefiles?
Must the fields be named the same for all 10 shapefiles?


texts: Google Map - saving latitude and longitude of route/path to my server to show in future the selected route/path and navigate accordingly

texts:
I am planning a software/mobile application where the tasks are like below.

An interface where a admin person will create a route map customisation by dragging a polygon path and also may add google suggested shortest path. So by combining both path software will save all latitude and longitude in server.
Later, when any user will go to that path he will see the customised path and can navigate.

We are using Google licensed version. My point is there is some privacy policy in Google like below

So my question is if I save latitude and longitude only does it violet the Google privacy policy?

If violet then what will be approaches to achieve my goal?


texts: Can ModelBuilder interface let user choose tool from among those within toolbox?

texts:
I would like the model user interface let the user select a tool (or submodel) from among a list.

Is it possible to prepare such a model using ModelBuilder without having to write a python script?


texts: Creating raster displaying graduated Distance to polygon or feature using QGIS

texts:
I am trying to create a raster which displays the relative distance of a pixel to a set of polygons of a specific value.

I have satellite derived fire scar data going back for the last 20 years.

For the purposes of my analysis any vegetation older that 3 years is classified as long un-burnt (this analysis is for a particularly fire prone region). The rest is classed as young vegetation.

It is easy enough to produce a shp/raster which, displays these two classes either in separate files or together. My first thought was to convert the files to geotif and do a proximity analysis in QGIS, I then tried a multi distance buffer but this didn't work properly either.

The blue polygons are old vegetation, the pink are younger. Ideally the raster will display a graduation scale of distance from the blue polygons (including the blue polygons which would display 0m/km).


texts:
Now I'm using OpenLayers to integrate OpenStreetMap into my web. Besides, I want to hide/remove all the labels on map (country, city, road, lake...).

How can I do that without establishing my own map server?


texts: Mapserver display empty map when i add more then 1 raster file in 1 table postgis

texts:
Here is a part of mapfile with layer and WMS:

  WEB
      METADATA
          wms_title "Diplom"
          wms_abstract "no matter"
          wms_onlinesource "http://localhost/cgi-bin/mapserv.exe?map=D:/Diplom/sample.map&"
          wms_enable_request "*"
          wms_srs "EPS    G:4326"
      END
  LAYER
  NAME rastertable
  STATUS       ON
  TYPE         RASTER
  DATA "PG:host=localhost port=5432 dbname=Diplom user=postgres password=**** schema=public table=rastertable"
  PROCESSING "SCALE=AUTO"
  PROCESSING "NODATA=9999"
  END

I load 2 tif files to postgis db use raster2pgsql -a -t 1000x1000 1.tif 2.tif schema.table and psql query.

Mapserver 3.0.4
Postgresql 9.4.0.
Logfile from mapserver:

  [Fri Mar 20 15:31:30 2015].251000 CGI Request 1 on process 3480
  [Fri Mar 20 15:31:30 2015].257000 msDrawMap(): rendering using outputformat named png (AGG/PNG).
  [Fri Mar 20 15:31:30 2015].257000 msDrawMap(): WMS/WFS set-up and query, 0.000s
  [Fri Mar 20 15:31:30 2015].260000 msDrawRasterLayerLow(rastertable): entering.
  [Fri Mar 20 15:31:31 2015].263000 msDrawMap(): Layer 3 (rastertable), 1.006s
  [Fri Mar 20 15:31:31 2015].263000 msDrawMap(): Drawing Label Cache, 0.000s
  [Fri Mar 20 15:31:31 2015].263000 msDrawMap() total time: 1.012s
  [Fri Mar 20 15:31:31 2015].377000 msSaveImage(stdout) total time: 0.114s
  [Fri Mar 20 15:31:31 2015].378000 mapserv request processing time (msLoadMap not incl.): 1.127s
  [Fri Mar 20 15:31:31 2015].378000 msFreeMap(): freeing map at 021ACF58.

When i trying to use http://localhost/cgi-bin/mapserv.exe?map=D:/Diplom/sample.map&layer=rastertable&mode=map in my browser he is display white screen, but if i use Db Manager in QGIS my raster files displayed normally.


texts: Is it possible to have a multi-modal network in QGIS?

texts:
I would like to design a multimodal network in QGIS and connect it with real-time traffic data (for road, rail and inland waterway). Is it possible to have a multimodal network in QGIS? And how can I get the real-time traffic data in Austria?

I already have a multimodal network in ArcGIS and I know that there is DTF files from NAVTEQ for roads. Is there any similar file formats for rail and inland waterway?


texts:
I was trying to create a Basetool which will collect points onMouseDown event. I was trying my tool to draw a line with the points collected. To collect points I used List. I am just putting the code here.

            IPoint screenPoint;
            screenPoint = new ESRI.ArcGIS.Geometry.Point();
            screenPoint.X = X;
            screenPoint.Y = Y;

            IScreenDisplay screenDisplay = (m_application.Document as IMxDocument).ActiveView.ScreenDisplay;
            IPoint mapPoint  = screenDisplay.DisplayTransformation.ToMapPoint(X, Y);
            point.Add(mapPoint); // Point is List of IPoints

Using 2 MOuseDown, I want to collect 2 points and pass the putCoords method the points to create a Line. Can anyone please help me ?


texts:
I have to filter the whole Malta dataset with Osmosis only considering cycleways, footways (including used nodes of those ways) and all nodes with the tag value tourism and amenity.

I write the following, but get argument 4 and 5 error no name.

osmosis --read-xml file="malta-latest.osm.bz2" --tag filter accept-ways highway="footway","cycleway" --used-way --tag filter accept-nodes amenity --tag filter accept-nodes tourism --tag filter reject-relations --used-node --write-xml file="malta-filter.osm.bz2"

Many thanks for the help in advance


texts: Time-enabled capabilty and attribute table lost when published as Tiled Mapping Service

texts:
I initially published a feature class with a massive number of polygons with large vertices on ArcGIS Online as Feature Service. This feature class has been enabled Time-aware capability. When opened and viewed in Map Viewer, it drew unacceptably slowly but good thing is time was retained.

I was directed by some one on the SE-GIS to publish the feature class to Tiled Mapping Service instead of Feature Service. I did that; the features were drawn pretty fast but Time was lost and no attribute table was available. Is this a bug or time/attribute table just not supported by Tiled Mapping Services?


texts: Forcing Feature Compare tool to iterate through multiple files in specific order?

texts:
I am working with multiple layers and I have been asked to find the percentage of change between them: GEOM, ATTR and ALL. I know that there is the Feature Compare tool and the Table Compare tool (I can't use the Detect Feature Changes tool as I am not using lines - only polygons.) and I've tried using the "Iterate Feature Classes" in order to iterate through each layer in my gdb, but it randomly selects what layers to compare against eachother.

For example if I have layers FC1 FC2 FC3 FC4, I want to compare FC1 to FC2. Then compare FC2 to FC3, FC3 to FC4 and so on. What the iteration tool does however is something like > FC1 to FC3, FC4 to FC1, FC2 to FC2, FC1 to FC4.

Is there a way to force Feature Compare tool to iterate through multiple files in a specific order?

I've had the thought of trying to create a script tool that lists the feature classes, but my coding skills are weak, and I'm not sure how to make this a functional tool.

I have access to both ArcGIS Desktop 10.5.1 and ArcGIS Pro.

And just  to be clear, I'm not married to these tools, they were just what was suggested to me. If there is a better way to do this comparison I'm open!

I need to process upwards of 200 layers eventually, I'm just currently working with 4 layers.


texts:
I have a pretty standard table with a geom column and a gist index on it.
I would like to have a fast pagination on this table as it contains millions of rows.
The query involves a ST_Intersects:
SELECT id, geom
FROM products
WHERE
    ST_Intersects(
        geom,
        ST_MakeEnvelope(%s, %s, %s, %s, 4326)
    )

However, there is no ordering with ST_Intersects, and the OFFSET clause adds a too large penalty.
Is there a way to achieve a fast and robust pagination with spatial geometries?


texts:
I am not a coder. However, I am familiar with a lot of the most basic concepts. My background is in the traditional GIS software/mapping industry but I have become a lot more interested in where the industry is going and I have realized I need to keep up. My issue is that the GIS community has become dominated by coding and web mapping, and this being the future, I realize I need to get good at this. They say the best way to teach yourself to code is a project, and I've got one, so here goes nothing.

I have undertaken a project that requires doing something like this. To do this, I need to learn how to use MapBox GL JS from the bottom up - from what I can tell, this means learning the basics for javascript essentially. I've already started doing that, using the IntelliJ IDE. What I don't understand is how to jump from the basics of the javascript language to using something like Mapbox GL JS. I'm not even sure what program/IDE (I think I'm on the right track with IDEA) I am supposed to edit something like this in, or how to make it work on our website beyond pasting that stuff into the HTML. I have all of the data and the style tiles made for the map I want to create, I just need to figure out how to code.

Now, I assume this probably represents a fairly big undertaking given how much of a novice I am, but can someone give me a roadmap for how to figure this out? I've got the time to do it. I also realize this is likely a stupid question to those who know what they are doing - hopefully some day I'm one of them!


texts: Adding iteration in model to join fields from several shapefiles into final one?

texts:
I have one grid (Atl1x1.shp) to which I need to add 9 Fields from 9 different shapefiles which are also grids (that is, 1 column from each shp, all with the same name). These shapefiles have the same number of records and the same Primary ID field as my original shp, so I used the JOIN FIELD tool where that PrimaryID would be both my Input Join field and my Output Join field, the original Atlantic1x1.shp would be the Input Table and the 9 shps would be the Join Tables.

Since I was in a hurry, I used the Join Field toolbox 9 different times, using the output of one Join Field as the Input Table of the following one. This way I ended up with a shapefile similar to the original Atlantic1x1.shp, with 9 extra columns in the attribute table coming from the 9 shps I used as Join Tables.

Since I will have to do this again in the future, I am working on a model to do it automatically. I do not know much about Model Builder, but I guess there must be a way to do this by adding one of the iterators.

In the pic you can see what I did so far.

Atl1x1.shp is the original grid

FishSpp1x1 is the directory where the 9 shps are stored.

My problem is that I don´t know how to tell the model to use Atlantic1x1 as the original Input Table for the Join Field, but then use the updated shp with the extra column as the Input Table in the second repetition, then that updated shp (with already 2 extra columns) as Input for the third repetition, etc.

Is there an easy way of doing this?


texts: How can I upload CSV file in QGIS Desktop App?

texts:
Please see link :
http://www.qgistutorials.com/en/docs/building_a_python_plugin.html

I am using the same process for uploading a CSV file in a folder but it gives me an error.

I am using this code:

def select_output_file(self):

    filename = QFileDialog.getSaveFileName(self.dlg, "Select output file ","", '*.txt')
    self.dlg.lineEdit.setText(filename)

It gives the error:

  Couldn't load plugin Importtool due to an error when calling its classFactory() method

  Traceback (most recent call last):

  File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 219, in startPlugin  plugins[packageName] = package.classFactory(iface)

  File "C:/Users/n/.qgis2/python/plugins\Importtool__init__.py", line 34, in classFactory  from .Import_tool import Importtool
    File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 478, in _import

  mod = _builtin_import(name, globals, locals, fromlist, level)

  File "C:/Users/n/.qgis2/python/plugins\Importtool\Import_tool.py", line 71
      self.dlg.lineEdit.clear()

  IndentationError: unexpected indent

Please see below image.


texts: v.kriging Addon not installed - GRASS GIS 7

texts:
I am using GRASS GIS 7. When I try to install "v.kriging" module using "g.extension" command, it is not installed. Though it says that the extension is successfully installed, I get the error "'v.kriging' is not recognized as an internal or external command, operable program or batch file."

I then tried to get the code from the following link:

http://wingrass.fsv.cvut.cz/grass70/addons/grass-7.0.1RC2/

But all folders in v.kriging.zip are empty.

Can you help me with installing the extension on Windows 8.1?


texts: I have a data set from a collard animal separated by days, I want to create an individual line from each day of data to know how the animal has moved in the forest. When I use the plugin points2path it only creates lots of line but I can not differentiate each day, I am using wrong the plugin? Does any body know a better way to do it?


texts: Batch Converting ASCII into Feature Classes using shapefile lists

texts:
Problem.  I have 413 ASCII files totalling 4 billion points, and my system cannot handle converting them to MULTIPOINT so that I can generate Terrains.  Therefore I would like to create 12+ subsets of the area, using the lists generated from a shapefile(s).

I've made my seperate Feature Classes listing the variable number of files that will make up each subset.
I can generate (and print) a list of files from the Feature Class using the SearchCursor function.
I'm running into problems using the generated list(s) to actually create the MULTIPOINT Feature Classes.

Here is my code:

# Import modules
import arcpy, os, sys, traceback
from arcpy import env

# generate a list of ascii files from the area-subset feature classes
env.workspace=r"F:\RTAMS\Petawawa\11-Library\LEGEO_Aug2012\LEGEO_Aug2012.gdb\LiDAR_ASCII"
fc=arcpy.ListFeatureClasses()

for fc in fc:
    cur=arcpy.SearchCursor (fc)

    for row in cur:
        fileName=row.TILENAME

        #this is the location of the ASCII files we want to import, so we change the workspace
        env.workspace=r"F:\RTAMS\Petawawa\11-Library\LEGEO_Aug2012\LiDAR\UTM18_WGS84\ASCII\FEATURES"

        # we should now have a list of files, which we will use to import into area-subset MULTIPOINT files
        # we set up some parameters for the ASCII3D to Feature Class tool

        try:
            arcpy.CheckOutExtension("3D")
            outWorkspace=r"F:\RTAMS\Petawawa\11-Library\LEGEO_Aug2012\LEGEO_Aug2012.gdb\LiDAR_ASCII_to_MULTIPOINT"
            inFormat="GENERATE"
            outFC=os.path.join(outWorkspace,fc+"_a")
            outType="MULTIPOINT"
            # find the ascii files listed in the search cursor
            txtList=arcpy.ListFiles(fileName)
            #Execute the ASCII3D_to_Feature_Class
            arcpy.ASCII3DToFeatureClass_3d(txtList,inFormat,outFC,outType)
            print outFC + " is complete."

        except arcpy.ExecuteError:
            print arcpy.GetMessages()

        except:
            # Get the traceback object
            tb = sys.exc_info()[2]
            tbinfo = traceback.format_tb(tb)[0]
            # Concatenate error information into message string
            pymsg = 'PYTHON ERRORS:\nTraceback info:\n{0}\nError Info:\n{1}'\
                .format(tbinfo, str(sys.exc_info()[1]))
            msgs = 'ArcPy ERRORS:\n {0}\n'.format(arcpy.GetMessages(2))
            # Return python error messages for script tool or Python Window
            arcpy.AddError(pymsg)
            arcpy.AddError(msgs)

Hoping you can help me out.


texts: Convert PLSS (Section, Township, Range, Meridian) to Latitude and Longitude using R package or Python

texts:
I have data in the form of Section, Township, Range and meridian. I have 300,000 data points from all over the USA.
Is there any way i can write R or python scripts to batch convert all the points to latitude and longitude?
What kind of equations they use to convert it into latitude and longitude?

I found this online tool which gives latitude and longitude of centroid and four corners of the section, but I want to do the same thing using R or Python.

http://www.earthpoint.us/TownshipsSearchByDescription.aspx

I tried to search online but couldn't find any information regarding this except online tools where we put information for each point manually.

first few lines of my data:

Township    Range1  Section  QQ            Footages
154 N       96 W      3     SWNE       1980 FNL  1980 FEL
154 N       96 W      3     NWNE       660 FNL  1980 FEL
136 N       83 W     22     NENW       330 FNL  2310 FWL
141 N       81 W     18     SESE       660 FSL  660 FEL
133 N       75 W     35     NWSW       1980 FSL  660 FWL
132 N       48 W     11     NWSW       2310 FSL  330 FWL
157 N       85 W     16     SESE       60 FSL  300 FEL
140 N       77 W     6      SWSW       330 FSL  330 FWL


texts: Why QGIS/PostGIS makes Sliver / small polygon when using erase (difference) tool unlike ArcGIS (QGIS vs ArcMap)?

texts:
When erasing b polygon from a, sliver(or small) polygons occured along the a polygon boundary in QGIS(3.4.5).
But ArcMap(10.3) does not make any slivers at all. I don't know which difference between two program exists.
I am struggling not to make any polygons and to know the reason. But, I cannot find any issues.

Both are not ring-self intersection polygons (st_isvalidreason)
Real Gaps are not shown.

I perhaps considered 'Tolerance issue'.. but I have no idea...

[


texts:
I currently have an app with two menu bars. One being the AxToolbarControl in a UserControl and one in a separate menu bar containing in-house stuff.

I currently have the issue where if I select something in the AxToolBarControl and then select something in the other menu-bar the item in the AxControl is still checked.  This obviously results in two commands being executed on the mousedown event handler.

I've looked at the IControl and IToolbar but it looks like .Checked is a ReadOnly.

Any ideas?


texts:
As far as I can tell, Google Maps does not release transit data or provide a transit API. Is there a one stop shop for transit Data:

Transit stops
Routes
Fare hierarchy and schedule
Ridership (optional)
Routing logic (transfer node, terminator, park n ride, etc)
Transit type (scheduled / call n ride)
Transit mode (Bus, shuttle bus, accordion bus, high speed subway, LRT, tram, commuter train, high speed train)
Transit organization make up (Transit is served by more than one org)


texts:
I have this kind of raster which has no pixel count in his attribute table.

I'd like to calculate the pixel count by class:

These are the raster's properties


texts: QGIS 3: equivalent tool/plugin of ArcMap's Spatial Join within a Distance

texts:
I would like to do a join by location between two point layers in QGIS 3 (3.4.10, Coruña). I would like to get the mean (or median) value of an attribute all the points of point layer 2 within a 5 meter radius of every point of point layer 1.

In ArcMap my approach would be:

Point layer 1 & 2: spatial join,  fill in the field map as "Mean" for the attribute, set join as "within a distance" and fill in 5 meters.

My approch in QGIS3 would be:

point layer 1: Add autoincremental field so i have a unique id for every feature
point layer 1: Buffer of 5 meter
buffered layer & point layer 2: Join attributes by location (Summary), fields to summarise -> the attribute field en summaries to calculate-> mean
Point layer 1 and summarised layer: Join attributes by field value on the base of the id field made in step one.

I wal through this process all the time, so I was wondering if there is a stand alone tool/plugin in QGIS3 for this?

Thank you!


texts: Joining a In_memory table with an SDE table

texts:
I have a in memory table that gets created, and want to join it to a table from SDE.
I have tried the following code.

Private Sub JoinTables(sdetable As ITable, inMemoryTable As ITable, joinField As String)

++ Create the MemoryRelationshipClass that defines what is to be joined
Dim pMemRelClassFact As IMemoryRelationshipClassFactory
Set pMemRelClassFact = New MemoryRelationshipClassFactory
Dim pRelClass As IRelationshipClass
Set pRelClass = pMemRelClassFact.Open("test", sdetable, "ID_NUMBER", inMemoryTable, "ID_NUMBER", "forward", "backward", esriRelCardinalityOneToMany)
' ++ Perform the join
Dim pRelQueryTableFact As IRelQueryTableFactory
Dim pRelQueryTab As ITable
Set pRelQueryTableFact = New RelQueryTableFactory
Set pRelQueryTab = pRelQueryTableFact.Open(pRelClass, True, Nothing, Nothing, "", True, True)

  ' ++ Print the fields   Dim pCursor As ICursor   Set pCursor = pRelQueryTab.Search(Nothing, True)
     Dim pField As IField
Dim pFields As IFields
Dim intI As Integer, intJ As Integer
Set pFields = pCursor.Fields
intI = pFields.FieldCount - 1
For intJ = 0 To intI
Set pField = pFields.Field(intJ)
    Debug.Print pField.Name   Next intJ
End Sub

but i get a type mismatch error on the line

  Set pRelClass = pMemRelClassFact.Open("test", sdetable, "ID_NUMBER", inMemoryTable, "ID_NUMBER", "forward", "backward", esriRelCardinalityOneToMany)

Is it not possible to join a in memory table to an sde table, or is there a better way of doing this?

Thanks


texts:
I would like to convert all mxd's in a directory and it's subdirectories to ArcGIS 10 mxd's. This is so that I can then run python to automate some updates within them.

Is this possible?

This is my second post and I was also wondering if I should post sample code?

=== UPDATE -- Here's the code (for more usefully) converting 10 files to 9.3

#Downgrades all ver 10 mxds in directory and subdirectories to version 9.3
#and saves them with _93 extension.

import arcpy, glob, os

path = os.getcwd()
for pathname, directories, filenames in os.walk(path):
    for filename in filenames:
        if filename.lower().endswith(".mxd"):
            mxd = arcpy.mapping.MapDocument(os.path.join(pathname, filename))
            filename, ext = os.path.splitext(os.path.join(pathname, filename))
            if mxd.dateSaved:
                    print mxd.dateSaved
                    mxd.saveACopy(filename + "_93" + ext, "9.3")
                    currentMxd = filename + "_93" + ext
                    mapDoc=arcpy.mapping.MapDocument(currentMxd)
                    print "Updated: " + filename
            else:
                    print filename + "is not version 10 therefore not altered"


texts: Is there a silent installation command without the desktop shortcuts?

texts:
I installed QGIS with Microsoft SCCM 2012 and I want to make a silent install without desktop shortcuts.

I managed to do a silent install with / S but I found no s command to delete desktop shortcuts.

I have Windows 7 64 client PCs

Thank you


texts: Toggle geojson labels using control with leaflet

texts:
I am attempting to create my first leaflet map.  I am using OpenStreetMap tiles and adding geojson features as overlay layers. I have been following tutorials with much success, until I started trying to add labels.  I would simply like to add static labels to several features, and be able to toggle them on/off with the layer control.  I have been using leaflet.label.  Is this possible?

I have tried the following to give me labels with mixed results:

function onEachFeature(feature, layer) {
    if (feature.properties && feature.properties.name) {
          <!-- shows label only on mouseover -->
          //layer.bindLabel(feature.properties.name, { 'noHide': true});

          <!-- creates static labels, but remain when layers are turned off -->
          //var label = new L.Label();
          //label.setContent(feature.properties.name);
          //label.setLatLng(layer.getBounds().getCenter());
          //map.showLabel(label);
       }
   }

Here is my layer control:

var overlayMaps = {
    "Layer 1": L.geoJson(layer1geojson, {onEachFeature: onEachFeature}).addTo(map),
    "Layer 2": L.geoJson(layer2geojson, {onEachFeature: onEachFeature}).addTo(map),
    "Buildings": L.geoJson(buildingsGeojson, {onEachFeature: onEachFeature}).addTo(map)
  };
var legend = L.control.layers(baseMaps, overlayMaps, {collapsed:false});
legend.addTo(map);

UPDATE

I used the following technique proposed by ghybs to add labels to geojson features and have the ability to turn them off with the layer control.

var buildingGroup = L.layerGroup().addTo(map);
L.geoJson(buildings, {
    style: buildingStyle,
    onEachFeature: function(feature, layer) {
        var label = new L.Label();
        label.setContent(feature.properties.name);
        label.setLatLng(layer.getBounds().getCenter());
        buildingGroup.addLayer(layer);
        buildingGroup.addLayer(label);
    }
});
// Repeat for any additional layers with labels

// Add Group to overlays
var overlays = {
    "Layer1": layer1Group,
    "Layer2": layer2Group,
    "Buildings": buildingGroup
};

// Add layer control to map
L.control.layers(baseMaps, overlays, {collapsed:false}).addTo(map);


texts:
Using ASTER DEM, I have created Flow Accumulation. However, to create stream order, I was asked to input flow accumultation and stream raster. But I cannot find how to get the stream raster. Some forums say that I should threshold flow accumulation. How do I do this? Please guide me in the right direction.

I am using ArcGIS 10.


texts: Multipart to Singlepart and Explode not working in ArcMap?

texts:
I am trying to separate multipart polygons, however in case of Multipart to Singlepart program seems to be working for around 15 minutes and then it crashes and when I try to Explode I get the error:"The Explode command could not be completed."

Is there any other way to separate such polygons?


texts: CRS defaults to WGS 84 on opening a project without OTF

texts:
In QGIS 2.18.15, I open a project and my component shapefiles have migrated to the west coast of Africa. Here's what's leading up to this:

I've disabled OTF on startup since there are well-documented geometry calculation issues that arise with OTF enabled.
My default projection is set to EPSG 102710 (NAD 1983 State Plane New Hampshire US Feet). All component shapefiles are projected to this as well.
When trying to calculate acreage for various polygons, these calculations are wonky, thus:

The "output preview" for the expression $area / 43560 correctly displays acreage, and yet:
The calculated column in the attribute table returns a value that is 1/10.7639 of the actual value. There are 10.7639 square feet in a square meter.

When I add a new vector layer, the "selected CRS" always displays as WGS 84 and I must manually select EPSG 102710, even though this is set to the default.

So it seems QGIS is defaulting to WGS 84 and metric measurements. If I turn OTF off, my shapefiles return from the African Atlantic to their rightful places in New Hampshire.

I do have several QuickmapServices layers included in my project.

Happy to provide other info if needed.


texts: Deleting features with date attribute before today's date using ArcGIS REST API?

texts:
I have a published feature service where road closure information is input to be shared with the public. The data input contains a couple of project date fields. I'd like to use a script that runs on a schedule on the server machine to automatically delete records before today's date. It looks like the way to properly do this is to use ESRI's ArcGIS REST API. We are running ArcGIS Server 10.1.

The date field is DISCONTINUEDATE.

I've worked out a query and it looks like I'd to delete all records WHERE:

DISCONTINUEDATE < SYSDATE

I'm brand new to the REST API and just not sure about the code to set up a delete using the above query.


texts: OpenLayers is not adding marker to my map

texts:
I have an OpenLayers map and I create it using a JavaScript function. I then have a button which when pressed calls another JavaScript function which has code to add a marker to the map.

For some reason the marker is not added.

Code for OpenLayers map creation:

function init()
    {
        var map, baseLayer;
        var options = {numZoomLevels: 5, isBaseLayer: true,};

        map = new OpenLayers.Map('map');
        baseLayer = new OpenLayers.Layer.Image("My map","map.png",
        new OpenLayers.Bounds(-2000, -2000, 1000, 1000),
        new OpenLayers.Size(1500,1500), options);
        map.addLayer(baseLayer);
        var center = new OpenLayers.LonLat(-100,-100);
        map.setCenter(center ,0);
        map.addControl(new OpenLayers.Control.LayerSwitcher());
        map.addControl(new OpenLayers.Control.MousePosition());
        var size = new OpenLayers.Size(21,25);
        var offset = new OpenLayers.Pixel(-(size.w/2), -size.h);
        var icon = new OpenLayers.Icon('http://www.openlayers.org/dev/img/marker.png',size,offset);
        var markers = new OpenLayers.Layer.Markers( "Marker" );
        map.addLayer(markers);
}

Code for adding marker:

function addMarker(x, y)
    {
        markers.addMarker(new OpenLayers.Marker(new OpenLayers.LonLat(x, y),icon));
    }

Any help is much appreciated :)


texts: QGIS creates a custom CRS when using a defined CRS

texts: I am trying to use the world_azimuthal_equidistant CRS in QGIS 3.4.2. When I try to re-project a layer from WGS 84 (the default) to the World_azimuthal_equidistant CRS, it automatically generates a custom CRS instead of the world_azimuthal_equidistant CRS. Is there a reason a custom CRS is automatically being generated?


texts: How to create polygon from my coastline in qgis?

texts:
I've downloaded the data from Natural Earth (coastline at a scale 50 million), I've cut the coastline into my study area, I wanted to clip other vector datasets as well, for this I wanted to convert my coastline to a polygon. I've tried 2 methods: 1) Polygonize and 2) Line to Polygon (both are in processing tool), but as you can see from both of the results of this method, it gave me weird polygon shape.

Are there any other tools to create my coastline to a polygon and how?

This is my study area, which is with 2 polygon

This is the result of "Polygonize" (coordinates: user-defined coordinates)

This is the result of "Line to polygon" (WGS84, 4326)


texts: Buffering a PointGeometry in arcpy only returns a line?

texts:
I am trying to buffer a PointGeometry object in arcpy to a Polygon geometry object with ArcGis 10.2.

The PointGeometry object is in lat/lon WGS84 coordinates and as far as I understand, I can only use PointGeometry.buffer by specifying the buffer in the same coordinates.
To specify the buffer in metres, I would need to use: arcpy.Buffer_analysis.

So I setup an empty geometry to do the buffer in like so:

spatialRef = u"GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],\
           PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];\
           -450359962737.05 -450359962737.05 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision"

pointGeom = arcpy.PointGeometry(arcpy.Point(-27.607361, 152.768311), spatialRef)
temPol = arcpy.Geometry()

pol = arcpy.Buffer_analysis(pointGeom, temPol, "50 meters")[0]

The resulting Polygon object however, is a long thin line, rather than the expected circular polygon.

Any ideas what I am doing wrong?


texts:
I've found ggmap's geocode() function really useful over the years as it's got a 'fuzzy' search implemented by Google. However, I'm wondering if there are alternatives, because

geocode() seems unreliable sometimes
overdependence on google
it doesn't use/support open source API solutions

I've seen the overpass  and geonames packages.

Has anyone else encountered this problem and discovered another solution, preferably using CRAN packages such as geonames?


texts: Adding custom icon doesn't work in Leaflet?

texts:
I'm trying to add custom icon to a map. All seems ok but when I click to add markers on my map I get the Leaflet icon by default and not the one I liked to have. Perhaps is in the way I get the IconUrl from a JSON or something else but I don't see where.

            displayMarkersByCat(offersFromCat,iconMarkers);
            console.log(offersFromCat);
            offersFromCat = [];
            displayIconByCAt(iconUrl);
            console.log(iconUrl);
            iconUrl = [];
        });
    });
}

function displayIconByCAt(iconUrl) {
    console.log(iconUrl);

    var data = iconUrl.toString(data);
    console.log(data);
    createNewicon(data);
}

function createNewicon(data){
    console.log(data);
    var iconMarkers =  new L.icon({
        iconUrl: data, // value of the URL img by cate
        iconSize: [38, 38], // size of the icon
        iconAnchor: [22, 94], // point of the icon which will correspond to marker's location
        popupAnchor: [-3, -76] // point from which the popup should open relative to the iconAncho
    })
    console.log({iconMarkers});
}

function displayMarkersByCat(offersFromCat,iconMarkers) {
    console.log(offersFromCat);
    console.log(data);
    // resetInterface();

    var markers =  new L.layerGroup();
    for (var i = 0; i < offersFromCat.length; i++) {
        var data = offersFromCat;

        console.log(data[i]);
        console.log(data[i][0]);
        marker = new L.marker([data[i][0], data[i][1]]).bindPopup(data[i][2]);

        markers.addLayer(marker,{icon : iconMarkers}).addTo(mymap);
        console.log(markers);
    }
    mymap.addLayer(markers);
}

And when I init my map for the first time I have static markers wich work well

var gastro = L.icon({
    iconUrl: 'www/images/alimentation.png',
    iconSize: [38, 38], // size of the icon

    iconAnchor: [22, 94], // point of the icon which will correspond to marker's location

    popupAnchor: [-3, -76] // point from which the popup should open relative to the iconAnchor

});
console.log(gastro);
var markerGastro = L.marker([43.6490449,
    0.5885573000000477], { icon: gastro }).addTo(mymap);

markerGastro.bindPopup("
I am a popup.").openPopup();

I just thinking it's where is the problem , my markers have been created already ????


texts: Plugin Reloader everytime i use custom plugin

texts:
Whenever I run my join csv script via custom plugin I have to run the plugin reloader for my script to take effect.

How can I make my script run without clicking plugin reloader everytime?

Heres my script:

from PyQt4.QtCore import *
import qgis
from qgis.core import QgsMapLayerRegistry, QgsVectorJoinInfo
import csv

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
  # Change to your shapefile name
  if layer.name() == "Naguilian":
    qgis.utils.iface.setActiveLayer(layer)
    shp = qgis.utils.iface.activeLayer()

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
  # Change to your csv name
  if layer.name() == "naguiliandb":
    qgis.utils.iface.setActiveLayer(layer)
    csv = qgis.utils.iface.activeLayer()

# Set up join parameters
shpField='PIN'
csvField='u_pin'
joinObject = QgsVectorJoinInfo()
joinObject.joinLayerId = csv.id()
joinObject.joinFieldName = csvField
joinObject.targetFieldName = shpField
joinObject.memoryCache = True
shp.addJoin(joinObject)


texts: How to smooth XY Co-ordinates of linear feature?

texts:
I need to smooth xy co-ordinates. see the image:

In this case, i need to remove 5,6,7,8 and 9th co-ordinates or need to make those co-ordinates  as a straight line between 4th and 5th. is there any logic to do this?


texts: Editing text in data view from layout view

texts:
I have text and arrows that I have drawn into my map in data view, using the drawing tool bar. I am often working in a fixed map scale, meaning when in layout view if I see text that needs to be adjusted I have to: go back into data view, pan over to the text, adjust it, and then go back to data view and pan back. I specifically DO NOT want to draw these items into layout for a wide range of reasons.

I can of course make my scale adjustable, never zoom in, and use the back arrow to readjust the view. There are lots of little tricks that help but if there is a tool that allows me to adjust the text in data view from layout view it'd save me a lot of time.

If this tool is part of some obscure extension that is okay I just want to know if it exists.

I am using ArcMap 10.5


texts:
NOAA provides QPF data on rainfall in the form of shapefiles. The shapefiles include .prj files with information on the coordinate system used:

GEOGCS["GCS_Sphere_EMEP",DATUM["D_Sphere_EMEP",SPHEROID["Sphere_EMEP",6371200.0,0.0]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

I was under the impression that NOAA data used WGS84, but it looks like this is not the case and it uses GCS_Sphere_EMEP. Am I misinterpreting the .prj file?

I am using this shapefile alongside route data from Google Maps which uses WGS84. I am not familiar with the Sphere_EMEP coordinate projection. Should I simply reproject the shapefile to use a WGS84 projection?


texts: Apportioning neighborhoods and census tracts - ArcMap

texts: I have a census tract shapefile and a neighborhood shapefile, and the end goal is to analyze vulnerable populations within each neighborhood using the data from the census tracts. Currently some of the census tracts overlap in to other neighborhoods, which I assume would lead to inaccurate results. I assume I need to apportion the neighborhood shapefile to the census tract in order to correctly analyze neighborhoods and particular vulnerable populations. I have seen a previous post on here titled Apportioning Neighborhood Information to Census Tracts using ArcGIS desktop?. A user described how to do it, but have been unsuccessful with these directions.


texts: Merge several .lyr files that contain a raster photo column

texts: Shapefiles do not support raster files, so I created a geodatabase and made feature classes of the points. I then added a column titled "photos" and added the photo as a raster to the point. I now have several .lyr files with rasters attached to the attribute table and ideally I would like to merge all of these layer files together so the photos (raster files) will all be contained in one .lyr file. I guess I should have merged the point shapefiles first, then made it into a feature class and add the raster files as photos to the merged file and created one .lyr file, but we didnt do it that way and I'm wondering if there is a shortcut to merge the .lyr files instead of attaching the raster files all over again. I am using Arc 10.3.


texts: QGis's Geoferencer plugin isn't able to paste text that has been added to clipboard via mac's "pbcopy" terminal command

texts:
I recently wrote a little script that used mac's pbcopy shell command to dump some text into my clipboard for use in the QGis program. For some reason I'm not able to paste stuff that has been pbcopy'd into my clipboard into the Geoferencer plugin on qgis. I recently discovered that I can paste the text into other elements of the native qgis program, but not the plugin. The text can similarly be pasted into any other program.

I think this has something to do with the way the Georeferencer plugin was written in Python, what with its somewhat unique at times deep bindings into the low levels of the operating system. Maybe some python devs could assist here.

I would save a lot of time in the long run if I could get the georeferencer program to respond to clipboard text that came from the pbcopy program.