texts: Seeking tools to calculate routing distance

texts:
I'm working to a web application, and I need to calculate the routing distance (not the distance "as the crow flies") between a starting address and several final addresses. I don't need to render the routes on the map.

The application will have a lot of queries, so I can't use Google APIs and probably the best solution is to perform the task not with an external service, but on a local server.

I think the only way is to use OpenStreetMap.

Which tool can I use to perform the calculation?

At the moment geocoding (700 conversions daily) is performed by Google Maps and Nokia's HERE APIs. Do you recommend me some other tools also for this?
The number of daily point-to-point routing calculations are approximately 100,000


texts:
I'm a heavy ArcGIS user, and was wondering how one may bring shapefiles onto CARTO.

I see that CARTO does not allow one to upload .shp datasets. Is there a easy conversion anyone would recommend?


texts: Make collector change symbology after registering an object

texts:
I use ArcMap 10.4.1 and the latest Arcgis Collector.

I have been using Arcgis Collector for some time now, but I still don't know how to make the symbology of an object change.

What I want, is for a point feature class that is already in the project (lets say houses) to change after it has been inspected. From a red square to a green square. I've been told that this is possible to do, but I can't find a single tutorial for this.


texts:
Based on other questions about how to "georeference" an image file, I took an SVG from our graphic designers, saved as a .dxf in Inkscape, and then imported into QGIS. I'm new to QGIS (we typically use ArcGIS but I can't import vectors there), so I could be trying to do something that another tool would be better for.

The image above was taken after numerous attempts. I import the .dxf, turn it into a GeoPackage, give it the same projection as the Esri basemap I'm using, turn on editing, and then use Vector Bender. I have tried using just 1 point translations, 2 point translations/scales/rotation, 3 point affines, and 3+ point bends. I have tried these separately and in different combinations. I cannot seem to figure out the right set of steps to get the black outline to actually overlay with the United States geographically.
I know you cannot see it well that far zoomed out, but the black outline is a more simplified/"graphical" representation of the states as well as counties within Ohio. My end goal here is to get that into a shapefile that I can use to more quickly plot points/distances when we're outputting data requests.


texts: Convert X & Y columns to a SHAPE column (using a method that would perform well in Portal)

texts:
I have a non-spatial table in an external system. The table has X & Y columns.

The table will be synced to an Oracle 18c/10.7.1 geodatabase table or FC (via the ArcGIS REST API). There will be roughly 200 inserts or updates from the external system per day.

I would like to serve up the synced geodatabase object (table or FC) to a Portal web map as a spatial point layer.

Question:

Is it possible to convert the X & Y columns to a SHAPE column -- using a method that would perform well in Portal?

Examples:

A feature class with an empty SDE.ST_GEOMETRY column. When rows are inserted or updated, a custom trigger would populate the value in the shape column (via st_point).
An XY event layer.
A registered spatial view that generates a calculated SDE.ST_GEOMETRY column (via the st_point function).
A registered spatial view that generates a calculated SDO_GEOMETRY column.
Something else?


texts: Searching a value from lists, get repeated results and append results to list in Google Earth Engine

texts:
Before, I asked that for search or filter values in two-segment list in Google Earth Engine at this link(Search or filter values in two-segment list in Google Earth Engine) and this question is answered by Daniel Wiell with very good explanation (https://gis.stackexchange.com/a/392774/112517).
I want to search for value in first column of the list, get finding result from second column and save to result to another new list. I want to search value from first column of lists and get second column value as result, for example, I want to search 101 and I get 12155 and save 12155 to new list.
var list = ee.List([
  [101, 12155],
  [123, 12333],
  [102, 12157],
  [102, 12222],
  [102, 12244],
  [103, 13333],
  [103, 12158],
  ]) // This list contains thousands of paired values

If this list have repeated values such(102), I don't know how to achieve other paired values(12222, 12244). For example, if I search for 102, I get only 12157 as result.
How can I achieve other paired values(12222, 12244) and append results to new list?


texts: Disable Leaflet WMS getFeatureInfo popup while using the Measure tool

texts:
I have a simple Leaflet map to which I've added two functionalities:

measure tool using leaflet.measure plugin;
WMS getFeatureInfo popup using the L.TileLayer.BetterWMS plugin;

The default map behaviour in this context shows both popups - measure popup and info popup.
Screenshot:

I have added the code to jsFiddle.

So my question is:

How can I prevent the info popup from showing up while using the measure tool?

Edit

I managed to disable the popup by using

$( ".leaflet-control-measure" ).click(function() {
  map.off('click', layer.getFeatureInfo, layer);
});

but now I need to turn them back on. I tried using jQuery toggle(), but something weird happens to the measure button.

$(".leaflet-control-measure").toggle(function() {
  map.off('click', layer.getFeatureInfo, layer);
 } , function() {
  map.on('click', layer.getFeatureInfo, layer);
});


texts: I need to load .shp files in Hive more precisely in HortonWorks sandbox through Beeswax, can anyone provide example of Hive commands or steps for Beeswax ?


texts: Loop using variable - r.mapcalc GRASS GIS

texts:
I am doing threshold classification of NDVI rasters. One equation looks like that:

r.mapcalc 'prog_r43_n65=if(ndvi_r43_n65>=0.6,1,2)'

Now I want to make a loop for more than 2000 rasters, so I created this:

for red in {43..64}
do
for nir in {65..193}
do
progi="ndvi_r"$red"_n"$nir    #I have files named ex. ndvi_r43_n65
progowa="prog_r"$red"_n"$nir
r.mapcalc '"$progowa"=if("$progi">=0.6,1,2)'
done
done

And here is the problem.. It does not read "progi" as a variable.. Does somebody guess why?


texts: Cannot center IActiveView at IPoint in C#.NET add in

texts:
I'm testing a simple portion of my add in and can't get it to center/pan to a specified point, I've included what is currently in my OnClick() method. No errors appear it just refreshes and nothing happens (10.1 SP1, C#.NET VS 2010).

IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
IActiveView actView = mxdoc.FocusMap as IActiveView;
IPageLayout pLayout = actView as IPageLayout;
ISpatialReferenceFactory srFactory = new SpatialReferenceEnvironmentClass();
int gcsCode = (int)esriSRGeoCSType.esriSRGeoCS_NAD1983;
ISpatialReference spatRef = srFactory.CreateGeographicCoordinateSystem(gcsCode);

try
{
    IPoint point = new PointClass();
    double xCoord = 33.217918;
    double yCoord = -98.158749;
    point.PutCoords(xCoord, yCoord);
    point.Project(spatRef);
    actView.Extent.CenterAt(point);
    actView.Refresh();
}
catch (Exception ex)
{
    MessageBox.Show(ex.Message);
}


texts:
I'm running QGis master v1.9.0 on xubuntu and the 'fetch python plugins' option is absent. This guy had this problem with vanilla 1.7.5 but fixed it by upgrading to 1.9.0 master. Nishant had this problem (version unknown); not sure if he fixed it but of the two solutions posted, the first (sudo apt-get install python-qgis) is already installed, and the second (uninstall, reboot, reinstall) changed nothing.

Any ideas?

Arguably irrelevant but actually not: is there somewhere we can find the estimated landing date for QGis 2.0? If it's due tomorrow, for example, there's little point trying to debug a problem in the 1.9 master only for all such problems to be eradicated with the fresh 2.0 install... :)

Thanks!!


texts: Formatting dependent parameters in Python Toolbox?

texts:
I am setting up my parameters for my tool in python toolbox. So far it looks like this:

The user inputs a tif into the Input Raster. I want the Select Bands (optional) tab to be a drop down menu that includes all of the bands contained in the .tif in the Input Raster. What changes do I need to make to my code (included below) to make this happen? I do not want to use the script wizard available in arcmap.  This has to be done in Python Toolbox only.

import arcpy,os,sys

class Toolbox(object):
def __init__(self):
    self.label =  "Exploding Rasters"
    self.alias  = "ER"

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

class ExplodeRasters(object):
    def __init__(self):
    self.label       = "Explode Rasters"
    self.description = " This tool will take an input multiband raster and extract its individual bands.  "

def getParameterInfo(self):

# Input Features parameter
in_features = arcpy.Parameter(
    displayName="Input Raster",
    name="in_features",
    datatype="DEFeatureClass",
    parameterType="Required",
    direction="Input")

# create select all bands button
all_bands = arcpy.Parameter(
    displayName="Select all bands",
    name="Select_All",
    datatype="Boolean",
    parameterType="Required",
    direction="Input")
TIFF = arcpy.Parameter(
    displayName="TIFF",
    name ="TIFF",
    datatype="Boolean",
    parameterType="Optional",
    direction="Input")
JPG = arcpy.Parameter(
    displayName="JPG",
    name="JPG",
    datatype="Boolean",
    parameterType="Optional",
    direction="Input")
IMG = arcpy.Parameter(
    displayName="IMG",
    name="IMG",
    datatype="Boolean",
    parameterType="Optional",
    direction="Input")

# Range of Desired Bands parameter
select_bands = arcpy.Parameter(
    displayName="Select Bands",
    name="Select_Bands",
    datatype="DEFeatureClass",
    parameterType="Optional",
    direction="Input")

select_bands.parameterDependencies = [in_features.name]
select_bands.schema.clone = True

#select bands for NDVI band parameters
NDVI_red= arcpy.Parameter(
    displayName="Select a red band for NDVI calculation",
    name="NDVI_bRed",
    datatype="DERasterDataset",
    parameterType="Optional",
    direction="Input")

NDVI_NIR= arcpy.Parameter(
    displayName="Select a NIR band for NDVI calculation",
    name="NDVI_BNIR",
    datatype="DERasterDataset",
    parameterType="Optional",
    direction="Input")

#out directory
OutDir=arcpy.Parameter(
    displayName="Output Workspace",
    name="Out_Directory",
    datatype="DEWorkspace",
    parameterType="Required",
    direction="Input")

# Create filename output prefix
prefix = arcpy.Parameter(
    displayName="Output filename prefix",
    name="output_prefix",
    datatype="String",
    parameterType="Required",
    direction="Input")

parameters = [in_features, all_bands, select_bands, NDVI_red, NDVI_NIR, prefix, OutDir, TIFF, JPG, IMG]

return parameters

def execute(self, parameters, messages):
    Out_Dir=parameters[6].valueAsText
    in_raster=parameters[0].valueAsText
    out_prefix=parameters[5].valueAsText

messages.addMessage("INPUT RASTER=" +in_raster)
messages.addMessage("Yay, you're doing great!")

#individual bands within multiband rasters already exist, so overwrite

arcpy.env.workspace=in_raster
arcpy.env.overwriteOutput = True

# get a list of the bands that make up the raster
bRng = arcpy.ListRasters()

# loop through the bands and export each one with CopyRaster
for ThisBnd in bRng:
    InBand  = ThisBnd
    bndDesc = arcpy.Describe(InBand)
    NoData  = bndDesc.noDataValue

    outRaster = os.path.join(Out_Dir, out_prefix + ThisBnd+".tif")

    #copy single band rasters
    arcpy.CopyRaster_management(InBand,outRaster, format = "TIFF", nodata_value = NoData)
return


texts: Preventing ArcGIS Online from truncating long field names?

texts:
According to the File Geodatabase Size Limits document, a field name can be up to 64 characters long:

  Field name length: 64 characters

I have a standalone fGDB table containing a column with a name 38 chars wide. The field's alias is the same as the field name.

I add this to an ArcMap document containing a single polygon layer, and publish the MXD to ArcGIS Online. The ArcGIS Online REST endpoint for the standalone table shows that the field name has been truncated to 31 characters long.

Is this a bug with ArcGIS Online, or is there a workaround?


texts: PyQGIS QgsRasterCalculator output format

texts:
I am currently scripting a simple raster operation. But, since I must do it thousands of time I would prefer if the output could be in the same format as the input or, at least in Int16.

Yet, it seems it is not the case and the output weight is far superior to what I would expect based on the size of the two initial rasters.

calc = QgsRasterCalculator( 'Int(rast1@1  - rast2@1)', output, 'GTiff', resp1.extent(), resp1.width(), resp1.height(), entries)

Both raster (rast1 and rast2) are integer (Int16), yet the result is in float32.

Is there a way to compress it directly in QgsRasterCalculator?


texts: Merging NDVI classified and unsupervised classified image?

texts:
I have clipped forest area from an image and apply NDVI classification on forest area and in the remaining, I have done unsupervised classification.

So, I want to know after NDVI and unsupervised classification, how do I mosaic or merge them to get a single classified image for further processing?


texts: Getting catalog path from IFeatureClass using ArcObjects?

texts:
VB.NET against ArcGIS 10.1
I have a function which searches a geodb for a feature class name. If found I would like to display where it was found. Currently I only return the FeatureClass object like so:

Dim fcTest As IFeatureClass = FindFeatureClassByName(pWorkspace, fcName)

It works great but now I would like to display the full catalog path of the feature class object. Is that possible?

I've been looking for hours but can't seem to do it. The Feature Class could exist in a Feature Dataset. So the feature class could be in locations like

    E:\Batch\Delivered.gdb\Bridges
    D:\Data\Final\Infrastructure.gdb\EastValley\powerlines
    C:/projects/RedRiverBasin/data.mdb/streams
    C:/projects/Airports/USA.mdb/West/lax

Does that info get included in the FeatureClass object or do I have to tweak my function?

I tried:

Dim pDataset As IDataset = CType(fcTest, IDataset)

But pDataset.Name is only the name of the feature class and not the full catalog path and name including any Feature Dataset if that's where its located.


texts:
I have a WMS layer in my leaflet map:

var wmsLayer = L.tileLayer.wms('http://localhost:8080/geoserver/volonteri/wms?', {
            layers: "trgovine",
            transparent: true,
            format: 'image/png'
        }).addTo(map);

In backend, there is a simple PHP script that stores new points into PostGIS table which is served via geoserver as WMS.

I use an AJAX request to send data to backend PHP script. I need to refresh the layer after successful insert to draw the new point on map, but cannot find any method for refreshing. The code looks like this:

$.ajax({
        type: 'POST',
        url: 'http://localhost:8000/api/trgovina',
        data: {
            'name': name,
            'shop': shop,
            'geom': geometry
        },
        success: function(){
            /* Redraw/refresh/reload should go here */
        }
    });

How can I solve this? I remember the older versions of openlayers had something simple like .redraw() method, but I'm pretty new to leaflet.


texts: Calculate cover loss between years using r.mapcalc

texts:
I am working on forest cover raster data for 2 years and am interested in areas with >5% loss in forest cover. I am pretty new in GRASS GIS, so any suggestions for how to create a difference map of two years using r.mapcalc?

by the way, I have used zonal statistics in qgis, to calculate % cover in each year.


texts: OpenLayers Vector Layers: adding and removing features VS. redraw/refresh

texts:
I am currently using addFeatures and removeAllFeatures in my code to get the desired results and this is working OK.

However, I have been asked to ensure that the vector layers are not being recreated each time. "They should be loading new data not creating new layers." I'm in no position to argue, and I really think the way I have it set up via removing/adding is pretty decent.

The type of data is airplane data in an AJAX call. The airplane moves along its flight path.

3 layers in question: the airplane layer, the projected flight path, and the travelled flight path.

Currently I have something like so:

if (vectors_airplane != "" || vectors_projected_path != "" || vectors_flown_path != "") {
  vectors_projected_path.removeAllFeatures([projected_path]);
  vectors_flown_path.removeAllFeatures([tail_path]);
  vectors_airplane.removeAllFeatures([airplane]);
}

and then for adding (within the same ajax call):

vectors_flown_path.addFeatures([tail_path]);
vectors_projected_path.addFeatures([projected_path]);
vectors_airplane.addFeatures([airplane]);

These are vector layers but there is no associated WFS. It seems to use the refresh strategy appropriately there must be a WFS Protocol (please correct me if I am wrong). I don't want to add a WFS protocol to each layer here, as it seems unnecessary.

The goal: how to change this logic to not use removeAllFeatures and addFeatures and instead use refresh/redraw?


texts: I just imported some shapefiles to a PostGIS database, where can I find the files in my drive?

texts: When I installed PostGIS, the default directory for all data is C:\Program Files\PostgreSQL\13\data but I can't find any folder that seems to be related to the database (named 'discover') that I made.


texts: Georeferencing a numpy tensor (spatio-temporal data)

texts:
I have a CSV file in which each row represents a record of an event with the following columns:

latitude, longitude, time

I want to build a numpy tensor with n x m x 8760 (hours in a year) dimensions where each cell will represent the sum of the events in the corresponding space and time.

So, in order to correlate the CSV events with the numpy tensor, I need to georeference the numpy tensor in a spatial reference system (ideally WGS 84).

Until now I have defined the numpy tensor as follows:

import numpy as np
from numpy import array
# define a spatial grid with time dimension (8760 hours)
c = np.zeros((8760, 100, 100)).astype(int) + 1

My first five lines from my CSV file are presented below:

Longitude, Latitude, Date, Time
-1.60202910794,53.7518017092,3/17/2017,0815
-1.53349564433,53.7944009244,1/14/2017,1330
-1.5624701361,53.7642562523,1/1/2017,0805
-1.5624701361,53.7642562523,1/1/2017,0805

So how can I georeference a numpy tensor with n x m x 8760 ?


texts: How to implement predefined attribute filters in an Openlayers web map?

texts:
I imagine this situation: I've got relatively big point layer and I want to split(filter) the points to predefined selections.
Predefined selection would be based on the attribute(for example Name:John, Type:Highway etc.)
*Users must be abble to turn the selections on/off *

Think that something like wfs filter should be usefull

Better said wfs filter is great and I only want to put its functionality to a checkbox, button, form or whatever what user will be able to set up(his own way)

new OpenLayers.Filter.Comparison({
                        type: OpenLayers.Filter.Comparison.EQUAL_TO,
                        property: "TYPE",
                        value: "highway"}),

Where can I find something about this theme? Or how to solve it?

Thanks for all the advices.


texts: Mosaic of rasters is big, black rectangle in QGIS?

texts:
I am brand new to QGIS and the world of maps. I have four GeoTiffs inside QGIS. I can merge them just fine using the Merge command, but when I try to use the Mosaic command instead I just end up with a black screen and a new layer which has a Min value of NaN and a really low Max number.

Can anyone tell me why and how to fix that, please? I see someone else had a similar problem in ArcGIS:
Why is mosaic of rasters big, black rectangle in ArcGIS Desktop?

I've since followed this guide for calculating statistics and then re-run Mosaic but nothing changed, it's still just blackness.

(The reason I'm using Mosaic and not Merge is because in a bit I'm going to need this command with a 5th image that overlays the other 4, and this is just a simple test.)


texts: Georeferencing a Picture having Resolution and a nearly correct BBOX

texts:
Is it possible to get the Boundingbox of a Picture when I know the Resolution, the Height/width in pixels, the sum of Tiles in each row/column and a Coordinate that is in the cornertile but not the corner !

To make it a bit clearer, take a look at that picture: The green things are given(one coordinate in the cornertiles and the resolution) and the red pointer indicates what I want to derive. Each tile is 256x256 pixels.

Does anyone have an idea ?

yours,
Matthias


texts: How to copy selected features from one sde to another using ArcObjects?

texts: Assuming I have two SDEs or other datasources with the same datamodel. One contains features the other one is empty. In ArcMap I have layers from both datasources in my table of contents. Now I would like to copy some features that are already selected in the map view from one sde to the other using C# and ArcObjects. What would be the best way to do this?


texts: Tutorial to create Shapefile based on Districts using ArcGIS for Desktop?

texts:
I have been tasked to create shapefile of Amritsar, India (will be incorporated into R programming language)

So far, I see this 2-part video tutorial, where person is assigning numbers and icons, and drawing in roads, but I don't see where he is separating a place into districts.

I've been searching endlessly, please point me in the right direction

Most likely I will use ArcGIS for Desktop via a 2 month free evaluation.


texts: Postgis identify points on the polygon border

texts:
I have points geom in gid,pt_geom and polygon geom in gid,pl_geom. Some of the points are inside polygon and some are outside polygon, others are on the polygon boundary.

Is there any way to identify those points on the polygon by using PostGIS functions (e.g. ST_contains() or ST_Distance())?

Maybe Select points.gid from points inner join polygon on st_contains(polygon.the_geom,points.the_geom).


texts: How to obtain common boundary lengths per year using Cshapes? (1945-2016)

texts:
How can I obtain border length per year between neighbor countries using Cshapes? (1945-2016)?
My current approach relies on QGIS' virtual layers, like the one below (source). However, after running the query, I obtain border lengths equal to zero in cases where they shouldn't. Am I missing something?
Cshapes: http://nils.weidmann.ws/projects/cshapes/shapefile.html
SELECT u11.*, u12.id || ', len:' || round(st_length(st_intersection(u11.geometry, u12.geometry)), 4) AS "neighbor_info"
FROM "unit_1" AS u11, "unit_1" AS u12
WHERE st_intersects(u11.geometry, u12.geometry) AND u11.id <> u12.id
ORDER BY u11.id ASC


texts: Update Leaflet circleMarker position with Android GPS

texts:
An Android code sends the latitude and longitude coords to a (html)Leaflet map (stored in Android Asset folder), which add a circleMarker to the actual poisiton.

function updateLocation(lat,lng,time){

map.panTo(new L.LatLng(lat,lng));

var circle = L.circleMarker([lat, lng], 10, {
    color: 'red',
    fillColor: '#f03',
    fillOpacity: 0.5
}).addTo(map);

I want to update the position, not to add to map again and again, when the GPS sends new positions. How can i fix it?


texts: Adding OpenWeatherMap WMS layer using PyQGIS?

texts:
I would like to add in QGIS real time weather data.

I found this link and this link and I tried to follow the example, but without any success.

I've added a WMS weather data layer for precipitation to the QGIS map, as follows:

First, the parameters for the service:

service = 'crs=EPSG:900913&dpiMode=7&featureCount=10&format=image/png&layers=precipitation&styles=&url=http://wms.openweathermap.org/service'

Next, create the raster layer, specifying WMS as the type:

rlayer = QgsRasterLayer(service, "precip", "wms")

Finally, add the precipitation layer to the map:

QgsMapLayerRegistry.instance().addMapLayers([rlayer])

the example didn't show any error in the code and no data was added to the weather layer either, I'm not sure why?


texts:
I am a GIS Programmer/Analyst in Vermont with formal training as an Environmental Scientist with GIS.  I have found throughout my GIS career the need to learn programming skills.  I currently work for a consulting firm where the majority of my tasks are automation using Python scripts and creating custom ArcGIS components using ArcObjects for VB.Net.  I am completely self taught in these technologies.

I very much enjoy the programming work I do and I would like to move toward becoming a professional GIS Software Developer.  Can anyone recommend an appropriate path to take toward this goal?  Specifically, I am interested in formal training to learn the best practices in Software Development that could be applied toward GIS software.


texts: How to use OR in a CON statement, in the ArcGIS Raster Calculator?

texts:
I'm trying to compare two rasters, and create an output based on a formula like:

if raster1 > 10 or raster2 < 20 then raster3 = 30

How can I specify the OR in the raster calculator's CON statement?

The pages Building expressions in Raster Calculator and Conditional evaluation with CON don't help, and "or" seems to be keyword that Google ignores....


texts: Account could not be used to sign in because it is not authorized to use ArcGIS Pro?

texts:
Yesterday afternoon (Friday), I renewed my Maintenance of ArcGIS for Desktop so that I could download and run ArcGIS 10.3 for Desktop and ArcGIS Pro 1.0.

Being in Australia, I was hoping that the order being placed by our local distributor would mean that this morning all entitlements would be unlocked.

Most of that seems to have worked because I was able to download and install ArcGIS 10.3 for Desktop and ArcGIS Pro 1.0 painlessly this morning.

I have also been able to start ArcMap 10.3 without a problem so the license I was using at 10.2 seems to have carried over nicely.

However, when I try to start ArcGIS Pro I am blocked by a message of:

  The account could not be used to sign in because it is not authorized
  to use ArcGIS Pro.

I have tried to trace what is astray by reading Technical Article 43739 which led me to a blog entitled Resolve recent problems with ArcGIS Pro authorizations, as well as the Manage Licenses page but when I (the Administrator) try to assign myself a license I get blocked by a Maximum of 0 Licensed Members in My Organization (see below). It may or may not be relevant but I participated in ArcGIS Pro Beta and skipped its Prerelease.

Does anyone know whether the provision of the necessary licensing to enable me to use ArcGIS Pro this weekend is something that must be worked through with the local distributor on Monday or should it have become available automatically once the downloads were unlocked by Esri?

This is a Commercial and not Home Use license so I think the answer to ArcGIS Pro 1 year home use says I'm not authorized will not apply.


texts: Getting Tkinter to work in PyCharm using OSGeo4W64's Python

texts:
I am creating an external QGIS program that builds virtual rasters using the processing algorithm within QGIS. I have this working, but it is not very convenient since you need to type in each file name. I would like to create a file dialog box using Tkinter.
In PyCharm the code is colored like the import is working, but then once I run the code I get an error stating that I do not have the module _tkinter. I look in the Python library in my OSGeo directory and I see a folder called tkinter.
I went to the OSGeo Shell and created a python environment in which I tired importing tkinter and got this...
>>> import tkinter
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "C:\OSGEO4~1\apps\Python37\lib\tkinter\__init__.py", line 36, in <module>
    import _tkinter # If this fails your Python may not be configured for Tk
ModuleNotFoundError: No module named '_tkinter'

Any ideas as to what is happening? Would installing pip and then installing tkinter through that work? I am not quite sure how to do that in OSGeo and fear messing something up so I have not tried.
If not, does anyone know of any alternative solution for my problem?


texts: Error adding qgis.org repository public key to apt keyring

texts:
The old QGIS key ring qgis-2017.gpg.key is expired and generates a warning during apt update today.

Following the QGIS instructions, I tried to update the key as follows:

wget -O - https://qgis.org/downloads/qgis-2019.gpg.key | gpg --import
gpg --fingerprint 51F523511C7028C3

But instead of the normal output, I get an error:

$ wget -O - https://qgis.org/downloads/qgis-2019.gpg.key | gpg --import
--2019-08-16 18:36:35--  https://qgis.org/downloads/qgis-2019.gpg.key
Resolving qgis.org (qgis.org)... 104.27.138.220, 104.27.139.220, 2606:4700:30::681b:8bdc, ...
Connecting to qgis.org (qgis.org)|104.27.138.220|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 1839 (1.8K) [application/pgp-keys]
Saving to: ‘STDOUT’

-                   100%[===================>]   1.80K  --.-KB/s    in 0s

2019-08-16 18:36:35 (33.3 MB/s) - written to stdout [1839/1839]

gpg: key 51F523511C7028C3: 1 signature not checked due to a missing key
gpg: no writable keyring found: Not found
gpg: error reading '[stdin]': General error
gpg: import from '[stdin]' failed: General error
gpg: Total number processed: 0

How can I fix the error please?


texts: Grass GIS "v.in.ogr": How to prevent crossing lines from beeing split where they don't share a common node?

texts:
I want to import a road network into grass gis to run network analyses. I use "v.in.ogr" to import my road shapefile. The problem is that Grass GIS seems to split every crossing line no matter if they share a common node or not.

For a road network this is problematic for bridges or tunnels.
How can I prevent the lines from being split? I'd like to preserve the line as in the following image:

Lines that belong together have the same id.

Which commands will lead to my desired result?


texts: How is distance measured in TIGER/Line data for roads?

texts:
So I'm working with some TIGER data for network analysis, and it comes with no speed or time attributes, so I was able to assign average speeds to the various road types, and I was going to use field calculator to populate the time field using the Speed and Shape_Length fields, but all of the lengths are very small, 0.00xxx. Is there something I need to do to convert them, or are they just measured in some strange unit of distance that I am not aware of?

This is what I have:


texts: QGIS: How to visualise according to amount of answers in one polygon?

texts:
I am working with QGIS. I managed to join a CSV file to a layer containing zip codes. I would now like to visualise the zip code areas according to how many survey answers I have from that one zip code area. I can't figure out how to do it. It did it once for me automatically, but I have no idea how.
The pictures demonstrates the shapefile of the zip codes and how many hits there might be from one zip code. I would like it to be darker where there are more answers and lighter where there are fewer.


texts: Output of RasterToPolygon not supported as shapefile

texts:
I made a python code for arcGis 9.3 which do the following step :

  1) Convert all the raster (.tif) in a folder to polygon shapefile in
  another folder

  2) Intersect all the converted shapefile polygon (previous) with a
  mask

The problem is that the intersect function doesn't work and it returns shapefile with no populate fields.
So I try manually to use the intersect Tool with the converted polygon create by my python code. And I have this error when trying to put the converted polygon shape in input :

PolygonName # Does not exist or is not supported

And when I Copy/paste this same polygon in another folder, then it works.

My code is as following :

import sys, string, os, arcgisscripting, glob
import arcgisscripting
gp = arcgisscripting.create()

# Set a default workspace
gp.workspace = "Z:\PPR\PPR76\PPRI_AUSTREBERTHE\Modelisation_Ruissellements\Tentative Macro Amelioration MNT"

# Set a toolbox
gp.toolbox = "management"

pathImages = gp.GetParameterAsText(0)
pathOutputConversion = gp.GetParameterAsText(1)
pathOutputIntersect = gp.GetParameterAsText(2)
masque = gp.GetParameterAsText(3)

# ON CONVERTI TOUTES LES IMAGES DU DOSSIER EN POLYGON

for infile in glob.glob(os.path.join(pathImages, '*.tif')):
    flnm = os.path.basename(infile)
    # Avec cette fonction : fileName contient le CHEMIN+nom du fichier sans l'estension et fileExtension contient l'extension (.qqch). Ici on voulait huste le nom du fichier sans le chemin donc on l'a applique a flnm
    fileName, fileExtension = os.path.splitext(flnm)

    # Local variables...
    Input_Raster= os.path.join(pathImages, infile)
    Output_Polygon_Features = pathOutputConversion + "\ " + fileName+ ".shp "

    # Process: RasterToPolygon...
    gp.RasterToPolygon_conversion(Input_Raster, Output_Polygon_Features, "SIMPLIFY", "")

gp.AddMessage(" FIN DE CONVERSION RASTER TO POLYGON ")

# ON INTERSECTE TOUS LES POLYGON IMAGES AVEC LE MASQUE DU BV

for infile in glob.glob(os.path.join(pathOutputConversion, '*.shp')):
  flnm = os.path.basename(infile)
  # Local variables...
  Output_Feature_Class = pathOutputIntersect + "\ " +flnm
  Input_Features = masque +" #;"+os.path.join(pathOutputConversion, infile)+" #"

  # Process: Intersect...
  gp.Intersect_analysis(Input_Features, Output_Feature_Class, "ALL", "", "INPUT")

Thanks for your help,


texts: QGIS - change of CRS raster - wrong scale

texts:
I have downloaded a raster file from Natural Earth Data. I would like to use it to make a map of Europe so I wanted to change the projection. I apply a new CRS to the layer and set the project CRS from the layer CRS. It allows me to do so but nothing really changes and the scale becomes completely wrong (a distance of 1000kms becomes 86m).

I am sure there is something basic I haven't understood (I am new to GIS and QGIS). Any suggestions?


texts: QGIS locking TAB files within the same organisation

texts:
I work in an office which uses Mapinfo, and all of our spatial datasets are saved as TAB files. A few of us also use QGIS to supplement our work. We're having issues where if someone has a TAB file open in QGIS, it locks the file so nobody else in the office has read access to it. It also breaks any FME workspace processes.

Any thoughts on how to get around this? We have a lot of files so we'd rather not duplicate everything as shapefiles.


texts: Cannot access unique value vector renderer for image in ArcMap?

texts: I have had to upgrade to ArcGIS 10.3, I tried to open a .tif image, but in the properties of the image I can not find the unique values option instead I found vector render?!! (look the picture attached). Can any one explain to me what is this for, and how can I find the unique values to my image to add legend to it and start processing.


texts: how the grass-gis extension r.diversity work?

texts:
I have found a nice paper about r.diversity extension. But I don't understand why... it returns me empty raster without errors.

I have tried to compute diversity on ndvi from spot6 image...

r.diversity --verbose --overwrite input=SPOT6_P_20151113@PERMANENT prefix=div_ size=3-11

a r.info on one of these readers give me

r.info map=div__pielou_size_3@unsup_class
 +----------------------------------------------------------------------------+
 | Map:      div__pielou_size_3@unsup_clas  Date: Wed Apr  5 17:55:41 2017    |
 | Mapset:   unsup_class                    Login of Creator: delaye          |
 | Location: spot                                                             |
 | DataBase: /home/delaye/Documents/futurSahel/grassDB                        |
 | Title:    div__pielou_size_3                                               |
 | Timestamp: none                                                            |
 |----------------------------------------------------------------------------|
 |                                                                            |
 |   Type of Map:  raster               Number of Categories: 0               |
 |   Data Type:    DCELL                                                      |
 |   Rows:         811                                                        |
 |   Columns:      825                                                        |
 |   Total Cells:  669075                                                     |
 |        Projection: UTM (zone 28)                                           |
 |            N: 1766650.21572733    S: 1765434.19215405   Res: 1.49941254    |
 |            E: 468559.15364472    W: 467321.45407192   Res: 1.50024191      |
 |   Range of data:    min = -nan  max = -nan                                 |
 |                                                                            |
 |   Data Description:                                                        |
 |    gnr par r.li.pielou                                                     |
 |                                                                            |
 |   Comments:                                                                |
 |    r.li.pielou --overwrite input="SPOT6_NDVI_20151113@unsup_class" conf\   |
 |    ig="conf_diversity_3" output="div__pielou_size_3"                       |
 |                                                                            |
 +----------------------------------------------------------------------------+

How it works ? Regarding the paper I don't need a conf file for r.li ...


texts:
I would like to join attributes by location using standalone PyQGIS but I'm getting error. I have QGIS 3.16 installed on my machine.
from qgis.core import *

QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.16/apps/qgis", True)
qgs = QgsApplication([], True)
QgsApplication.initQgis()

import processing
from processing.core.Processing import Processing
Processing.initialize()

processing.run("qgis:joinattributesbylocation", {
            'INPUT':'C:/Users/ja/Inne/Desktop/bdot10k_bubd_a/bdot10k_bubd_a.gpkg',
            'JOIN':'C:/Users/ja/Inne/Desktop/POLSKA/PRG_PunktyAdresowe_POLSKA.shp',
            'PREDICATE':[1,5],
            'JOIN_FIELDS':['PNA','SIMC_id','SIMC_nazwa','ULIC_nazwa','Numer'],
            'METHOD':1,
            'DISCARD_NONMATCHING':True,
            'PREFIX':'',
            'OUTPUT':'C:/Users/ja/Inne/Desktop/projekt2.gpkg'})

And I'm getting the error:

Traceback (most recent call last):   File "", line 1, in
   File
"C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\tools\general.py",
line 108, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)   File
"C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python/plugins\processing\core\Processing.py",
line 141, in runAlgorithm
raise QgsProcessingException(msg)
_core.QgsProcessingException: Error: Algorithm qgis:joinattributesbylocation not found

What is the problem?


texts: Labeling coordinate value at 4 corners of graticule in ArcGIS Desktop?

texts:
Is it possible to label automatically lat/long values at 4 corners: upper left, upper right, lower left, lower right of the graticule in Arcgis (10)?

I have not found the way to do that so far.


texts: How to get layer description using ArcObjects?

texts:
Assuming we have a ILayer object, I'm looking for how to get the layer description.

I mean the description text of the layer properties:

I didn't find it into the ArcObjects documentation.


texts: "Join Attributes by location (summary)" only with more than 50% of overlapping area

texts:
I have been using the "Join attributes by location (summary)" tool in QGIS to help create heatmaps. I am wondering if there is a way to only join a feature if it overlaps more than 50% inside the respective area. For example, the picture below:
Orange represents a count of 1 and red a count of 0. I would prefer that only the bottom polygon is orange because it contains the majority of the green feature.
Any ideas?


texts:
Does anyone know of a way to find the standard deviation of an attribute lying within any given polygon?

The plugin dialog "Join Attributes by Location" gives only the options Mean, Min, Max, and Sum.

I'm using QGIS 1.6.0.


texts: Convert X,Y coord to lane id and road coords

texts: Is there any way or function that can convert x,y coordinates to get s coordinate, road id and lane id? I need this info for making OpenScenario file, so want to convert x,y coord present in xodr.


texts: How to remove strange gdaldem hillshade artifacts?

texts:
I've downloaded 9 chunks of 1 arc-second NED in IMG format from the National Map.  I'm trying to create a hillshade from them.

I'm using:

GDAL 1.10.0, released 2013/04/24

First, I combine them and change the projection:

gdalwarp -t_srs EPSG:900913 n40w111/imgn40w111_1.img n40w112/imgn40w112_1.img \
n40w113/imgn40w113_1.img n41w111/imgn41w111_1.img n41w112/imgn41w112_1.img \
n41w113/imgn41w113_1.img n42w111/imgn42w111_1.img n42w112/imgn42w112_1.img \
n42w113/imgn42w113_1.img uinta-projected.tif

Then I create a hillshade from uinta-projected.tif:

gdaldem hillshade -compute_edges -co compress=lzw uinta-projected.tif uinta-hillshade.tif

However, when I take a look at my shiny new hillshade, it looks like this:

Does anybody have any ideas what might be causing this grid artifact?  I've tried using gdal_merge.py instead of gdalwarp and I end up with the same result.

gdalinfo uinta-projected.tif:

Driver: GTiff/GeoTIFF
Files: uinta-projected.tif
Size is 9253, 12173
Coordinate System is:
PROJCS["Google Maps Global Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-12579287.992128284648061,5161229.105774102732539)
Pixel Size = (36.130094040215752,-36.130094040215752)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-12579287.992, 5161229.106) (113d 0' 6.00"W, 42d11'34.83"N)
Lower Left  (-12579287.992, 4721417.471) (113d 0' 6.00"W, 39d11'11.36"N)
Upper Right (-12244976.232, 5161229.106) (109d59'54.57"W, 42d11'34.83"N)
Lower Right (-12244976.232, 4721417.471) (109d59'54.57"W, 39d11'11.36"N)
Center      (-12412132.112, 4941323.288) (111d30' 0.29"W, 40d42'24.63"N)
Band 1 Block=9253x1 Type=Float32, ColorInterp=Gray


texts:
Like many people in geology and archaeology I need to work with local grids that have eastings and northings from a specific reference point, and typically have a difference alignment/direction of north to true north.

1)  Any way to create such grids in QGIS?

2)  If I need to rotate layers that contain groups of points, lines or polygons in QGIS is the affine translation plug in the only way? If so can someone point me to an online calculator that can covert simple rotation angle, scale factor and rotation orgin values into the transformation matrix values required by the plug in ([x'=, x+, y+][y'=, x+, y+])?

3)  Questions about local grids and rotation in QGIS are very common. How much time would it take a programmer to put together a plugin to create local grids and/or do vector rotations on groups of points, lines or polygons?

If answering – please note I am not a mathematician or programmer, and I have no idea what Python is.


texts:
I am using Django 1.3, PostGIS 2.0 , Postgresql 9.0 and GeoDjango.

I had a database with a table for places. This table had latitude and longitude for places in the US. I enabled this database with PostGIS. I added PointField to my model. Then, updated this table with

SELECT AddGeometryColumn('table', 'poly', 4326, 'POINT', 2)

CREATE INDEX "table_poly_id" ON "table" USING GIST ( "poly" )

UPDATE table
SET poly = ST_PointFromText('POINT(' || longitude ||' '|| latitude ||')', 4326)

Now, I want to make query to find places within 5 miles of given latitude and longitude of a user's location.
So, I did this

ref_pnt = fromstr("POINT(%s %s)" % ("-87.627778", "41.881944"))
records = Table.objects.filter(poly__distance_lte = (ref_pnt,D(mi=5)))

But, I am getting an error

DatabaseError: Coordinate values are out of range [-180 -90, 180 90] for GEOGRAPHY type

I also get an error when I do this

 SELECT ST_PointFromText(poly) FROM table

ERROR:  parse error - invalid geometry
HINT:  "0101000020E61000007E" <-- parse error at position 20 within geometry
CONTEXT:  SQL function "st_pointfromtext" statement 1

What could be the problem?


texts: Can you do aggregate functions among nearby geometry in PostGIS?

texts:
I want to do aggregation based not on group but on geometry, is that possible in PostGIS?

For example to calculate max, min, average price for all listings that's nearby a geometry point, how can I do something like:

select max(price), min(price), avg(price)
from price_table
where (ST_DWITHIN(geom_point, ST_SETSRID(ST_MAKEPOINT(-119.697321, 34.40592), 4326), 0.72))

If PostGIS doesn't have this feature, what other system should I look into?


texts: Using numeric variable in where clause of ArcPy SearchCursor?

texts:
I know this question has been answered but I am not seeing an adaptation that I can understand. I am trying to use a numeric variable in a where clause and just having a hard time figuring it out.

fc1 = "L:\ENGINEERING\collind\ArcGIS\Default.gdb\sedaSampleHistory"

fieldname = "pocid"

delimfield = arcpy.AddFieldDelimiters(fc1, fieldname)

pocidCnt = 3

cursor = arcpy.da.SearchCursor(fc1, ["pocid", "cl2total"], delimfield + " = pocidCnt ")

I have tried all kinds of ways to get pocidCnt to work in the where clause and I am not having any luck. If I write it as

cursor = arcpy.da.SearchCursor(fc1, ["pocid", "cl2total"], delimfield + " = 3 ")

it works fine.

What am I doing wrong?


texts: Why am I getting duplicate ghost points (features), i.e., one point in the feature class but two points?

texts:
When editing some of the attributes in a point feature class in ArcMap 9.3, duplicate labels started showing up.  When I click on the point, I discovered that the point has duplicated itself BUT there is only one line in the feature class.  If I highlight the row (the feature), I find that there are two points in one row - almost like a ghost point.  I can't get rid of them and it seems that it gets worse (more duplicate ghost points) the more I edit.  If I select all in the attribute table, the number selected is the correct number but the table adds the ghost points to the total number of records.  Confussed?  So am I.  Tried stop editing and deleting version and restarting -- starts back at normal and then starts replicating ghost points as I edit.

Thoughts????
Thanks,
Sandy


texts: Define a circle extent using the circle tool - arcgis javascript api

texts:
I  have the following script for drawing a circle and I want to limit the user to a predefined radius:

$('#circle').click(function () {
       clearPreviousSearch();
       map.disableMapNavigation();
       $('#run').fadeOut();
       buttonclickvalue = 'circle';
       drawToolbar.activate(esri.toolbars.Draw.CIRCLE);
         });

Where should I specify that I want to define the specific radius using the ArcGIS API script?:

var symbol = new SimpleFillSymbol().setColor(null).outline.setColor("blue");
       var gl = new GraphicsLayer({ id: "circles" });
       var geodesic = dom.byId("geodesic");
       map.addLayer(gl);
       map.on("click", function(e) {
         var radius = map.extent.getWidth() / 100;
         var circle = new Circle({
           center: e.mapPoint,
           geodesic: domAttr.get(geodesic, "checked"),
           radius: radius
           });


texts: How to get India OSM map offline in geoserver which looks similar to online OSM map?

texts:
I need to OSM map of India having same look and feel as of Online OSM map. I have done following steps for the same :

I have downloaded India map OSM data from below link :

http://download.geofabrik.de/asia/india.html - >india-latest.osm.pbf
Exported this osm pbf file data in postgis DB using osmosis. so, below are the tables generated in DB :

actions, nodes, relation_members,  relations, schema_info, users, way_nodes, ways
I tried to publish these tables as vector layer in geoserver. however, it doesnt display anything.

So, please suggest me the next step/correct me if I am doing anything wrong.

My final goal is to get offline map of India (without internet),  which looks same as of online OSM map.

Based on other posts I tried using tilemill/mapbox. however, I am not clear about the imported India Map data in Postgis+tilemill/mapbox relation. And, the relation between tilemill/mapbox + geoserver.

I am implementing GIS solution using Geoserver+PostGIS+Openlayers technologies.


texts: ArcGis layer query: Average annual precipitation for area in U.S

texts:
I need to get average precipitation for certain land parcels. I've found several layers which could offer this information but queries do not return results or indicate error. Below is example of my request ( our application uses javascript api). I need result as number inches/yard or similar?

Query: Average annual precipitation (inches/yr) (ID: 21)

query that returns error
Query: Annual Precipitation (ID: 0)


texts: I imported the shapefiles from a government source but it looks like the land masses arent't aligned properly even though the crs on all layers are the same.


texts: Extracting part of string from field in QGIS Field Calculator?

texts:
I want to create a new field in my attribute table. My goal: returning all digits between the / and _ characters of the field named "location" by applying the next expression:
regexp_substr( "Text", '/(\\d*)\_' )
based on this working example: Obtaining specific part of string from field in QGIS attribute table?
However, the expression does not give the desired result in my case:

Does anyone know any way to modify/accommodate the expression?


texts:
I have a RasterLayer and a SpatialPolygonDataFrame.

I used the mask function to create a new RasterLayer.

The output is correct but it does include lot's of NA values.

Now I want to get rid of the NA values.

Is there an efficient way to solve this problem?


texts: How to get points within boundaries which overlaps the anti meridian

texts:
Within my application I have a map and I'm showing markers on the map which represent objects that I get from the database. The objects have latitude and longitude. The objects are selected only if their latitude and longitude are within the boundaries of the map. To be precise, it's a Google Maps block with a map and on every drag/zoom I request the server for a list of new objects passing the 4 values of the map boundaries that I receive from the Google Map. In GMaps JS API it looks like this:

bounds = map.getBounds(); // Bounds of the map.
NE = bounds.getNorthEast(); // North East corner.
SW = bounds.getSouthWest(); // South West corner.

After that, on every event fired, I send a request to the server sending this data:

{
    'x_min': SW.lng(),
    'y_min': SW.lat(),
    'x_max': NE.lng(),
    'y_max': NE.lat(),
}

It is Django/Python on the server, so I handle the request like so:

polygon = Polygon.from_bbox((x_min, y_min, x_max, y_max))
queryset = queryset.filter(geo_location__contained=polygon)

Now, all of that works fine unless... the right side of the boundary is crossing the anti meridian. So, the left side is for example 90 and the right side is for example -170. In this case I obviously get no results on the map.

I'm pretty much sure that's a common problem, but can't find any viable solution.


texts:
I know Line dash not working OL3 (Openlayers 3) has been asked earlier but the solution proposed is not working.

My code is :

var calculateStyle= new ol.style.Style({
      fill: new ol.style.Fill({
    color: [255, 255, 255, 0.6],
    opacity: 0.3
  }),
  stroke: new ol.style.Stroke({
    color: [0, 153, 255, 1],
    linedash: [40,40],
    width: 3
  })
});

I tried both [4,4] and [40,40] suggested in the question linked. But it's not working.


texts: LayerToKML output doesn't render on Google Earth

texts:
For the last few weeks I've been making .kmz files using LayerToKML_Conversion with no problem. They open in Google Earth perfectly.

Today, they suddenly stopped rendering at all. When I click the .kmz file, Google Earth opens and renders the map, but my KML doesn't appear. There are no error messages. The name of the .kmz file appears in 'Temporary Places' like normal, but doesn't actually appear (it should zoom in to Tanzania and render some yellow outlines):

I haven't changed anything at all in my code, it's just two lines of text:

import arcpy
arcpy.LayerToKML_conversion('PreHex280_merged', 'C:/Users/alell/Desktop/PreHAN_Control_hex280/final.kmz')

It started rendering when I added is_composite=True to the code, but it renders in a very low-res fashion that is not very useful.


texts: Simple homemade route optimization algorithm

texts:
I'm trying to do a simple route optimization without any fancy frameworks and insert-science-professor-name-here algorithms.
I have a startDepot and 100k waypoints of which I want to find about 25 waypoints, which are reasonable to visit in a single route. The accuracy doesn't have to be perfect.
Currently I'm searching for the nearest point to the start location, then search the nearest to the first point, then nearest to the previously found point etc, always excluding the previously found points. Now, the algorithm gives quite good results and even produces a reasonable route, but I really would like to include an endDepot also. Sometimes it may match the startDepot, sometimes not.
Any ideas of how to find the next point that prefers the direction to the endDepot?

I'm using postgis but the question is more related to algorithm generation


texts: My @media query breaking the state of the map in ArcGIS Javascrpit API

texts:
Anybody knows why my map is not displaying again after resizing my window?  I am using the next media query to hide the map when the width is less than 767px.

@media (max-width: 767px) {
.map { display: none !important; }
}

The map is hiding but is not restoring again when I maximize the window. If I use { opacity: 0} it's working. Why display:none is creating an inconsistent state of the map?

I need help, thanks

PS: sorry for my English.


texts: All Python scripts broken in ArcMap 10.7.1

texts:
I recently was attempting to run a Python script that requires the lxml package. I couldn't figure out how to use pip to install it as it always defaults to Python 3.x. So, as a work around I just copied the folders Python made when it installed it for Python 3 and just pasted them into the same folder for Python 2.7 for ArcMap. Well, now all the Python scripts just show a red circle with an X and when I run anything it doesn't run. It just pops up an error message with no text, just a yellow triangle with a ! in it.
Is there a fix or do I need to just reinstall?


texts: How to properly import QImage from PyQt4.QtGui?

texts:
I am starting with QGIS python console. Naturally, I came across the PyQGIS Developer Cookbook. While trying to cut my way through it, I encountered a lot of problems with "undefined attributes and names".

Some of them are caused by not importing some library. I had the same problem as described here.

  NameError: name 'QImage' is not defined

Even though I imported the PyQt4.QtGui library. At least, I thought:

Python 2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)] on N-PC
## Type help(iface) for more info and list of methods.
help(QImage)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
NameError: name 'QImage' is not defined
import PyQt4.QtGui
help(QImage)
Traceback (most recent call last):
  File "<input>", line 1, in <module>
NameError: name 'QImage' is not defined
from PyQt4.QtGui import *
help(QImage)
Help on class QImage in module PyQt4.QtGui:

class QImage(QPaintDevice)
 |  QImage()
 |  QImage(QSize, QImage.Format)
 |  QImage(int, int, QImage.Format)

And the help goes on. Why did QGIS fail to load QImage by the 1st import command? Or what actually happened here? I am using QGIS 2.4.0, and I get the same issue in 2.0.1.


texts: How to find points inside each polygon in Mysql?

texts:
I have points table (with latitude and longitude columns) and polygons table with geometry stored in Mysql. I want to get points inside one polygon, for example in polygon which has column name="California".

This is my query which returns all points inside all polygons, but I need all points inside polygon with name 'California':

$sql = "SELECT points.name FROM polygons, points WHERE ST_CONTAINS(polygons.geom, Point(points.longitude, points.latitude));"


texts:
I am having hard times trying to figure out how to possibly convert a dataset they gave me in a custom SR into a well known EPSG 32632.

The situation is:

I received a IFC whose coordinates were unknown to me.
I imported it into ArcGIS Pro using Data Interoperability letting the tool "guess" the SR for me
In my scene, the dataset (now a multipatch featureclass) is placed where it should be (in Vicenza, Italy) although Z coordinate is not reliable (but as of now I don't mind). I can confirm this by overlapping the dataset to a Esri Basemap
I read the SR is recognized as "IFC COORDSYS 0"
If I set the map SR to EPSG:32632, everything surprisingly stay in place and I don't have the Z issue anymore

Now, I need this dataset to actually be in 32632, but, no matter what I try, the result is all the time a strange "spear-like shape" while it should be a tunnel, or "disappear" from the mapo in 32632.

Here is what I have been trying so far:

using define projection to overwrite the prj of my multipatch from "IFC COORDSYS 0" to "WGS_1984_UTM_Zone_32N"
creating an empty multipatch featureclass from scartch assigning the 32632 SR and then copying and pasting the geometries from the original FC to this brand new one
importing the IFC specifying the 32632 as my SR

I then decided to save the custom projection to a prj and inspecting it to hopefully have a clue of what's going on.

Below are the two prj: the one that I am trying to project from ("IFC COORDSYS 0")*****, and the one I am trying to project to ("WGS_1984_UTM_Zone_32N").

Can you give me a hint to possibly find out a solution to this problem?

***** prj is made through right-clicking the custom SR in the Map properties and then "Save as Projection File..."

IFC COORDSYS 0

PROJCS["IFC_COORDSYS_0",
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["World_Geodetic_System_of_1984_GEM_10C",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Azimuthal_Equidistant"],
PARAMETER["false_easting",0.0],
PARAMETER["false_northing",0.0],
PARAMETER["central_meridian",-71.25807189916667],
PARAMETER["latitude_of_origin",42.41486358638889],
UNIT["MILLIMETER",0.001]]

WGS_1984_UTM_Zone_32N

PROJCS["WGS_1984_UTM_Zone_32N",
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",9.0],
PARAMETER["Scale_Factor",0.9996],
PARAMETER["Latitude_Of_Origin",0.0],
UNIT["Meter",1.0],
AUTHORITY["EPSG",32632]]


texts: Platforms for repetitive automated mapping?

texts:
I could give a lot of examples, but imagine you have to do, for example, a study in 10 different locations and have to generate a lot of maps with the same layout, legend.

Maps can contain vector and raster data.

It should be scriptable, because a large number of maps have to be made/remade in the future.

Which platforms can do this?


texts: Why many in-car navigation systems choose to go through the city center?

texts:
Most of the in-car navigation systems I've been using prefer to select an itinerary that goes through the city center, rather than getting around. And this is the case while the "Fastest" (not "Shortest") route option is selected.

Of course, mathematically, the shortest path from A to B is a straight line. So even though a city is between A and B, the shortest path is going through the center of that city, distance wise.

But (besides traffic, that requires to get real-time data, and is not the point of this post) cities have usually traffic lights, lower speed limit etc... that increase significantly the time to travel from A to B.

So the question, that I've been wondering for some time, why do many in-car systems go through city centers

Because even today's algorithms do not take into account the number of traffic lights to cross, city speed limits etc...?
Because navigation systems have been "required" to get people do shopping etc... within cities?
other reasons...


texts: Shortest distance between polygon and polyline

texts:
I have a shapefile with several hundred polygons (archaeological sites) and am trying to calculate the shortest distance for each to a single polyline (modern coastline). As part of this process, I am also trying to produce a new shapefile with polylines describing this shortest distance that I can then turn into annotated measured distance lines during map production.

I would prefer to use QGIS, but do have access to multiple other GIS platforms.


texts:
The code below is not working: This code doesn't affect the table at all

layer = point_layer3
selected_feature = layer.selectedFeatures()
print(len(selected_feature))
layer.startEditing()

for feature in selected_feature:
    a = f['parkeer_approx']

    b = parking_capacity[f.id()-int(1)]['parking_capacity']
    feature["parkeer_approx"] = a - b
    layer.updateFeature(feature)

layer.commitChanges()


texts: QGIS Modeler: Extract by location, multiple predicates

texts:
I'm creating a processing model and would like to use the 'extract by location' tool,

The problem is that when I choose this tool from the 'Algorithms' in the model creator, I cant find the double criteria as we can find in the one in the QGIS toolbox (normal use). The model creator gives an 'extract by location' with only one criterion.

Is there any way to have many criteria in the same query?


texts: Line thicknesses for PDF have no difference in ArcMap

texts:
I have a map with several layers with different thickness. When I export into PDF there is no difference in thickness. Practically lines 0.25mm and 0.5mm are exported all as 0.5mm. When I export this map into JPG the thickness it is OK, it can be visible even without zooming in the map.
The fast solution is, I can export into .svg and create the PDF by using Inkscape.
I need to export directly from ArcMap 10.5 because I want to automate some stuff.
This screenshot is of two images exported from the same layout.


texts: How to assign symbology in shapefile A based on data in shapefile B - QGIS 3.0.2

texts:
I have two shapefiles: "States" - contains polygons for the 50 US states, and "Widget_stores" - contains points where widget stores are located.  I want to adjust the symbology of "States" to change colors based on the number of stores in each state (ex. 0 stores, 1-5 stores, 6-10 stores, 10+ stores). Both shapefiles contain a STCODE field with the two-letter abbreviation for the state (ex. CA, FL, TX, MN, etc.)

I tried using the Field Calculator in the "States" layer to look in the "Widget_stores" layer and count the number of points (stores) located in each state either by geometry or by STCODE and assign that number to a new integer field NUMSTORES in the "States" layer so I can classify each state for the symbology.  I constructed an expression using the "aggregate" function because the documentation says it returns an aggregate value from another layer, however the result was not what I expected - for each of the 50 states it returned the total number of stores.  My expression must be missing a filter or some other operation needed to parse state-by-state.    Here's my expression:

  aggregate('Widget_stores','count_distinct', "STCODE")

Looking for suggestions on how to construct the proper expression!


texts: Spatial join with polygon and multipolygon mixed feature class

texts:
I have a feature class which contains 300,000 polygon and multipolygon features in the same fc. The goal is to get a result layer of spatial join (must to be one to many join) between this feature class to itself and during the process I have to only use open source GIS software, better code than geoprocessing tools. I tried:

GeoPandas sjoin - doesn't have one to many option.
QGIS 1 - Tried to read the fc in QGIS then use join by location. The join crushed after 20%.
QGIS 2 - Tried to convert fc to shapefile then use multipart to singlepart than dissolve to get the first mix of polygon and multipolygon. Again the join crushed in the same time.

Mabye QGIS cannot handle this kind of geometries mix?


texts: Calculating population density per sq km in QGIS Field Calculator?

texts:
I'm using QGIS 2.10 and I have a table with the following data:

Where "ALAND10" is the area in square meters and "DP0010001" is the total population.

How to calculate population density per 1sq km? I have ~30 000 entries?


texts: Positional argument error using rasterio to reduce raster pixel size

texts:
I'm new with rasterio and I need to reduce the pixel size of a GeoTiff. Using the code below
def resample(input_raster: str, output_raster: str = None):
    upscale_factor = 2
    with rasterio.open(input_raster) as src:

        # resample data to target shape
        data = src.read(
            out_shape=(
                src.count,
                int(src.height * upscale_factor),
                int(src.width * upscale_factor)
            ),
            resampling=Resampling.bilinear
        )

        # scale image transform
        transform = src.transform * src.transform.scale(
            (src.width / data.shape[-1]),
            (src.height / data.shape[-2])
        )

        with rasterio.open(
            output_raster,
            "w",
            driver='GTiff',
            height=data.shape[2],
            width=data.shape[1],
            count=1,
            dtype=src.dtypes[0],
            crs=src.crs,
            transform=transform,
        ) as dest:
            print(dest)
            dest.write(1, data)

I see this error:

test_raster.py:28 (test_resample) tmp_path =
PosixPath('/tmp/pytest-of-maxdragonheart/pytest-0/test_resample0')
def test_resample(tmp_path: Path) -> None:
    output_data = str(tmp_path / "resampled.tif")

   resample(
        input_raster=input_raster,
        output_raster=output_data
    )

test_raster.py:32:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  ../../raster.py:123: in resample
dest.write(1, data)

??? E   rasterio.errors.InvalidArrayError: Positional argument arr must be an array-like object
rasterio/_io.pyx:1338: InvalidArrayError

How I can solve?


texts:
I have just done a fresh install of ArcGIS Desktop 10.4.1 on my PC (new PC, never had ArcGIS on it before).  I tried to run the Union tool and it gave me an error

  ERROR 001143: Background server threw an exception.
  ERROR 001143: Background server threw an exception.

(yes it said it twice)

Re-running the tool gave the same error (several attempts).

The help doesn't really give much info as to what may have caused it.

Running with background geoprocessing disabled as suggested in the help is successful, but when enabled again it fails.

I hadn't yet installed the 64bit Background Geoprocessing, so I've installed that and now the tool works fine without error.

What are the requirements to run Background Geoprocessing in ArcMap?  Does Background Geoprocessing in ArcGIS require the 64bit Geoprocessing to be installed to run from within ArcMap, or is this just coincidence that it now runs (perhaps due to the restart of ArcMap?


texts: Combine aerial imagery pictures to create a layer

texts:
We are looking for software or a web service that we can put in pictures we took using our weather balloon to create a single image, that we can use as a layer for our GIS project. Preferably free, or free for education.

EDIT: We use a ESRI system


texts: Join leads to null value outside of edit session

texts:
I am using Arcmap 10.2 on Windows 7

I joined a shapefile with an Excel file, and the attributes are shown as  until I start an edit session. Then, the values are as they should be. They match up by verifying manually the attribute table and the Excel file. I clicked "Save edit" and then ended the edit session, at which point the attributes went back to .

Why is it doing that? How can I fix it?


texts: What is the max value for node ids in overpass?

texts:
I wanted to use overpass with custom data using the .osm format.

Overpass does not consider negative node IDs, so I added a prefix, staying <63 bits.

ex. : node id = "540000000000001" (in the .osm file)

But I have this error when I want to import data :

Reading XML file ... finished reading nodes. Flushing to database
..terminate called after throwing an instance of 'std::bad_alloc'
what():  std::bad_alloc
./bin/init_osm3s.sh: line 44:  2007 Broken pipe
bunzip2 < $PLANET_FILE 2008 Aborted
(core dumped) | $EXEC_DIR/bin/update_database --db-dir=$DB_DIR/$META $COMPRESSION

If I delete a few zeros, it works (ex. : node id = 9900000000001).

I've read that de value of node ids must be encoded in 64bits (signed). It seems that overpass allows less. What is its max value ?


texts: OpenLayers ol.proj.get returns null after adding projection definition to Proj4

texts:
I just stumbled upon the following problem :

OpenLayers is unable to retrieve a projection definition even though it has been previously registered with Proj4JS. Code below :

   import proj4 from 'proj4/dist/proj4'; // Proj4JS version 2.4.4

   import proj from 'ol/proj'; // OpenLayers version 4.6.5

   proj4.defs('EPSG:2154', '+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3 +x_0=700000 +y_0=6600000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs');

   console.log(proj.get('EPSG:2154')); // displays null

I can not figure out anything but a version mismatch as reported in OpenLayers 3 not recognising proj4 custom projection


texts: Exporting from Excel to csv and import in QGIS retaining encoding/special characters (with workaround)?

texts:
Problem:
When exporting Turkish character encoding into a csv file from an Excel Table (Excel 2013)  as a source file I can´t import the resulting csv into QGIS to still displaying the characters correctly .
I´m using: "Add Layer" - "Add Delimited Text Layer" and it doesn´t work no matter whether I use as encoding option "UTF-8" or "System".

Result is always:

This is possibly rather an excel problem since the problem concerns exporting from Excel. If anyone has any suggestion that would be great, otherwise anyone else reading this with the same problem could try the following workaround.

One possible workaround I use is using Access as an in-between-step:

Importing the Excel-Table into Access retains the special characters.
Exporting from Access into a txt or csv file you have the special button called "Advanced options" in the lower left of the export window where you can specify the Encoding of the resulting exported file. If I use "UTF-8" here and then import the resulting csv or txt file into QGIS selecting "Encoding "UTF-8" everything works fine.

Another Question on this site concerns exporting the other way round: Exporting from QGIS to Excel retaining special characters, see here:
Export shapefile with "non standard" characters to csv and read in excel

So how do I get Excel to export the characters correctly for QGIS to read them?

My System is
Windows 7 Professional 64-bit (6.1, Build 7601)
Service Pack 1
Language: German (Regional Setting: German)

QGIS Version is 2.12.0 since I haven´t had time to upgrade.


texts: How to rotate multiple polygons about their respective anchor points

texts:
Possible Duplicate:
  Rotate a feature in ArcGIS 10

I have a set of polygons, which are to be rotated in eight directions, about their respective fixed anchor points.

I am able to do so for a single polygon using Affine (Scale, Rotate, Translation), but unable to do so for multiple polygons in the same layer.

Can anyone provide a python script, for the same.

Regards,
KSV


texts: WMS bounded overlay with OpenLayers and GeoServer

texts:
I'm trying to design a heatmap overlay bounded within four given coordinates with OpenLayers using WMS protocol to communicate with my GeoServer instance, where my data and SLD heatmap description reside.
I can basically do it but I incur in two problems I'd like to fix:

I have to set "singleTile: false" in the WMS options, and this causes the heatmap to be recomputed for each sent tile, and thus have some discontinuity.
Here is an example of what I mean with discontinuity
I select bounds using maxExtent option, which causes entire tiles to be displayed, even if I'm only interested in a small portion of the selected tile.

Here is my code from the javascript side:
  format = 'image/png';
  heatmap = new OpenLayers.Layer.WMS(
"Rainfall", "http://localhost:8080/geoserver/postgis/wms",
{
  LAYERS: 'postgis:rainfall',
  format: format,
  transparent: true
},
{
  singleTile: false,
  visibility: true,
  buffer: 0,
  maxExtent: new OpenLayers.Bounds(-15000000, -6000000, -8000000, 2000000),
  isBaseLayer: false
}

My questions are:
What is the right way to do what I'm trying to do?
Is it possible to use WMS protocol to do that?, should I do it at a lower level(maybe in the SLD)?


texts:
We have a very complicated setup for a project we are doing, and it involves 40+ developers and two oracle databases, one with sde and one only with oracle spatial and it's quite big.

We'll here are my issues:

In the SDE database, ArcSDE is creating my linear tables as Linestrings, while in the other Oracle database, we have Multilinestring tables. I'm not sure how that got by, and there is a ton of code built on top of the oracle only schema.

How can I transform a simple linestring in a multilinestring? Here is what we've tried:

Using oracle append function to see if it returned a multi* geometric field. it does return a multi* geometry, but some parameters in the metadata are different
using sdelayer -o delete and sdelayer -o register to force a multilinestring. The procedure was fine, but even after new inserts, the oracle db is still seeing linestrings.
Merging two geometries inside arcgis. Still linesting geometry.

Any thoughts? Hacks? Solutions?

Changing the datamodel or code, specially code inside our oracle only db is an issue. Not possible.

Thanks for your attention and insights!

george


texts: How to determine latitude longitude point from a pixel within a MODIS tile?

texts: I need to know latitude / longitude of each pixel in a MODIS tile image in HDF in order to reproject from sinusoidal (MODIS) to mercator (Google Maps). I Can read values from HDF in my C++ program, I can reproject with proj4 having latitude / longitude of a point. I found a MODLAND Tile calculator  with source code in c++ but it have too many dependences in Linux libraries and I'm programming a Windows aplication. Any suggestion will be appreciated


texts: Can I export a layer's field value to metadata ?

texts:
Is there a way around for me to export a layer that has the field name "Quadname" and export it to metadata to capture all the quads names? I am looking for that to work on my metadata using Arcgis 9.3.1.

Thanks


texts: Make the nodata area of a resampled orthophoto overview white?

texts:
I have a set of orthophoto images which are each fully filled with data, but when placed together the coverage does not form a rectangle so there are nodata areas.

In MapGuide, I need to be able to display these nodata areas as white so that my users don't have to waste black ink when printing.

At large scale, where I display original full resolution images directly, this is no problem.  I just set my map's background colour to white, and the areas where there are no orthophotos display the background.

For performance, I need to be able to merge all of these source images into a resampled composite overview image for display at smaller scales where more of the orthophotos are viewable at once.

I've been attempting to use GDAL to merge and resample the overview, but by default it appears to create the  resampled composite GeoTIFF tile with black in the nodata areas, and MapGuide does not allow me to set black as transparent on colour rasters.

Is there a way for me to efficiently get what I want?

I have provided the answer I ended up with using GDAL, but would love to see solutions using other image processing utilities and GIS applications, both open source and proprietary.


texts: Is there a QGIS plugin to drop pins on a repeating pattern in a raster?

texts:
I'm looking for a QGIS plugin that drops pins on a repeating pattern in an imperfect raster grid, for example a satellite image. Ideally, one would select a polygon to drop within, a raster, give it the rough direction of rows and an approximation rough row and column spacing, and it will automatically drop pins on a repeating pattern, accounting for some degree of uneven spacing. Is there any plugin that does this or something similar?
We're specifically trying to drop pins on individual plants on a farm.


texts: Passing list as reclassify input table parameter in qgis.reclassifybytable algorithm

texts:
Following up on Save and load the classes definition table in the QGIS reclassify by table tool, I decided to run the qgis.reclassifybytable algorithm using PyQGIS to reclassify an aspect raster.
Here's the code:
import numpy as np
import processing

l = np.hstack([np.arange(1,5,1),1])
t = np.hstack([np.arange(-90,360,90)[1:]])
t = np.repeat(t,2)
t = np.hstack([0,t,360])
L = list(np.insert(t,np.arange(2,len(t)+2,2),l))

input_layer = QgsProject.instance().mapLayersByName('Aspect')[0]
output_layer = '/tmp/reclassify_output.tif'
params = {
  'DATA_TYPE': 5,
  'INPUT_RASTER': input_layer,
  'NODATA_FOR_MISSING': False,
  'NO_DATA': -9999,
  'OUTPUT': output_layer,
  'RANGE_BOUNDARIES': 2,
  'RASTER_BAND': 1,
  'TABLE': L # <---- I'd like to pass a list here!
}

alg_id = 'qgis:reclassifybytable'
processing.run(alg_id, params)

Which gives this traceback:
Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 108, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 159, in runAlgorithm
    ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "/usr/share/qgis/python/plugins/processing/gui/AlgorithmExecutor.py", line 71, in execute
    results, ok = alg.run(parameters, context, feedback, {}, False)
_core.QgsProcessingException: Invalid output value:

I particularly do like the Invalid output value: with nothing after! And the OUTPUT value is actually perfectly valid.
It took me some time to actually figure out a way that works, but with which I'm not 100% satisfied; copy-pasting the value of the list L that is printed on the console when I simply barely hit "L" in the Python console to show the content of this variable and to naively put it where L is used in the params variable in my code:
>>> L
[0, 0, 1, 0, 90, 2, 90, 180, 3, 180, 270, 4, 270, 360, 1]

And this works.
So I thought that maybe, for some reason, the algorithm is waiting for a string at this position, so I also try to enter str(L) instead of simply L but then:
Traceback (most recent call last):
  File "/usr/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 24, in <module>
  File "/usr/share/qgis/python/plugins/processing/tools/general.py", line 108, in run
    return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
  File "/usr/share/qgis/python/plugins/processing/core/Processing.py", line 159, in runAlgorithm
    ret, results = execute(alg, parameters, context, feedback, catch_exceptions=False)
  File "/usr/share/qgis/python/plugins/processing/gui/AlgorithmExecutor.py", line 71, in execute
    results, ok = alg.run(parameters, context, feedback, {}, False)
_core.QgsProcessingException: Invalid value for minimum: [0.0

Sounds very strange. It looks like it stops after reading the first value of the list.
Does anybody know what's going on here and how could I programmatically pass my list to the function params instead of entering it manually?


texts: Creating GeoTiff from Tiff+WorldFile+Proj

texts:
I'm currently using the global mapper export function to convert my tiffs+worldfile+projection file to geotiff.  Is there a commandline tool (for Linux) that lets me do the same thing?
"gdaltranslate" doesn't seem to use the projection file.


texts: Rule of Thumb for setting smoothing parameter

texts: I'm trying to model some binary data (disease infections) using indicator kriging in the Geostatistical Wizard in Arc 10.3. For presentation purposes I'm using the smoothing option. Is there a rule of thumb or an equation to use to set the best smoothing factor? The default is 0.2, but is that necessarily the best? I've looked many places and have yet to find any actual guidelines this way.


texts: Converting coordinates to EPSG:4326 for Google Maps

texts:
I have a MultiLinestring, example coordinates, (698745.50201423 4210011.04002000) and I would like to convert them to something like (24.362 42.458) to use in Google Maps. When I try to get the SRID from the shapefile, it tells me it is an unknown (0) SRID.

I use SELECT UpdateGeometrySRID('roads','geom',4326); to set SRID in 4326 but again the format of coordinates is not changed.
After I tried:

ALTER TABLE roads
    ALTER COLUMN geom TYPE geometry(MultiLineString,4326) USING ST_Transform(geom,4326);

and

  ALTER TABLE roads
  ALTER COLUMN geom TYPE geometry(MULTILINESTRING, 4326)
    USING ST_SetSRID(geom,4326);

but again the format of the coordinates is still like (698745.50201423 4210011.04002000).


texts: Google ai-platform model definition for earthengine usage

texts:
I want to host a DNN-model on the Google AI-platform and use this model in the Google Earth for prediction.
However, I can not seem to configure the input-output definition properly. There is an example here, but without any further explanation.
from tensorflow.python.tools import saved_model_utils

meta_graph_def = saved_model_utils.get_meta_graph_def(config.MODEL_DIR, 'serve')
inputs = meta_graph_def.signature_def['serving_default'].inputs
outputs = meta_graph_def.signature_def['serving_default'].outputs

# Just get the first thing(s) from the serving signature def.  i.e. this
# model only has a single input and a single output.
input_name = None
for k,v in inputs.items():
  input_name = v.name
  break

output_name = None
for k,v in outputs.items():
  output_name = v.name
  break

# Make a dictionary that maps Earth Engine outputs and inputs to
# AI Platform inputs and outputs, respectively.
import json
input_dict = "'" + json.dumps({input_name: "array"}) + "'"
output_dict = "'" + json.dumps({output_name: "impervious"}) + "'"

# You need to set the project before using the model prepare command.
!earthengine set_project {PROJECT}
!earthengine model prepare --source_dir {config.MODEL_DIR} --dest_dir {config.EEIFIED_DIR} --input {input_dict} --output {output_dict}

Is there any mapping/explanation available how the input_dict relates to the access in the Earth Engine?
var model = ee.Model.fromAiPlatformPredictor({
    projectName: PROJECT,
    modelName: MODEL_NAME,
    version: VERSION_NAME,
    inputTileSize: [pixel_size, pixel_size],
    proj: ee.Projection('EPSG:4326').atScale(SCALE),
    fixInputProj: true,
    outputBands: {'land_classes': {
        type: ee.PixelType.float(),
        dimensions: 1
      }
    }
});

I fail to understand the mapping and can therefore not use a model of my own.


texts: ST_DWithin exponentially slow. Cannot find what I am doing wrong

texts:
PostGIS version: 3.1
PostgreSQL version: 12.3
The machine I am working with has: 126G RAM, 48 CPU cores

Info:
I am getting started with PostGIS.
My goal is to get all the matching data between two points.
lv.geopoint and sub.geopoint both are GEOGRAPHY Points (SRID: 4326) and have GIST indexes on them.
My sub SELECT returns about 3k lines, my ‘valeurs_foncieres’ table however has 14 000 000 lines.
I do have BTREE indexes on valeurs_foncieres.id, caracteristiques_2018.id, caracteristiques_2018.num_acc, usagers_2018.id, usagers_2018.num_acc, vehicules_2018.id, vehicules_2018.num_acc.
The problem:
The query gets exponentially slow as I increase the distance of ST_DWithin.

Precision 100: 2sec
Precision 1 000: 10sec
Precision 10 000: 6min

Here is the query:
SELECT
    DISTINCT(sub.num_acc),
    sub.geopoint,
    sub.id
FROM
    (
    SELECT
        DISTINCT(u.num_acc) AS unumacc, c.*
    FROM
        usagers_2018 u
    INNER JOIN vehicules_2018 v ON
        u.num_acc = v.num_acc
    INNER JOIN caracteristiques_2018 c ON
        u.num_acc = c.num_acc
    WHERE
        u.grav = '2'
    ORDER BY
        c.id
) AS sub
INNER JOIN valeurs_foncieres vf ON
    ST_DWithin(vf.geopoint,
    sub.geog,
    1000,
    FALSE);

Here is the EXPLAIN:
HashAggregate  (cost=265577998.10..265578004.81 rows=671 width=49)
  Group Key: c.num_acc, c.geopoint, c.id
  ->  Nested Loop  (cost=9948.38..264845621.97 rows=97650150 width=49)
        ->  Unique  (cost=9947.84..10316.67 rows=6706 width=170)
              ->  Sort  (cost=9947.84..9964.60 rows=6706 width=170)
                    Sort Key: c.id, u.num_acc, c.an, c.mois, c.jour, c.hrmn, c.lum, c.agg, c."int", c.atm, c.col, c.com, c.adr, c.gps, c.lat, c.long, c.dep, c.lat_gps, c.long_gps, c.geopoint, c.geog
                    ->  Gather  (cost=3200.48..9521.63 rows=6706 width=170)
                          Workers Planned: 1
                          ->  Nested Loop  (cost=2200.48..7851.03 rows=3945 width=170)
                                Join Filter: ((u.num_acc)::text = (v.num_acc)::text)
                                ->  Parallel Hash Join  (cost=2200.06..6686.70 rows=2075 width=170)
                                      Hash Cond: ((c.num_acc)::text = (u.num_acc)::text)
                                      ->  Parallel Seq Scan on caracteristiques_2018 c  (cost=0.00..2859.90 rows=33990 width=157)
                                      ->  Parallel Hash  (cost=2174.12..2174.12 rows=2075 width=13)
                                            ->  Parallel Seq Scan on usagers_2018 u  (cost=0.00..2174.12 rows=2075 width=13)
                                                  Filter: ((grav)::text = '2'::text)
                                ->  Index Only Scan using vehicules_2018_num_acc_idx on vehicules_2018 v  (cost=0.42..0.54 rows=2 width=13)
                                      Index Cond: (num_acc = (c.num_acc)::text)
        ->  Index Scan using valeurs_foncieres_geopoint_idx on valeurs_foncieres vf  (cost=0.54..39477.72 rows=1456 width=32)
              Index Cond: (geopoint && _st_expand(c.geog, '1000'::double precision))
              Filter: st_dwithin(geopoint, c.geog, '1000'::double precision, false)
JIT:
  Functions: 30
  Options: Inlining true, Optimization true, Expressions true, Deforming true

Questions:
Is this normal? How can I decrease the execution time?


texts: Accessing Python dictionary by field value in QgsExpression

texts:
I'm writing a QGIS processing plugin. In this plugin, I use the processing algorithm pointsalonglines, and for the distance parameter, I'd like to use an expression which gives back the value of a dictionary with the key equal to a field value. Not sure that is clear but basically:
mydic = {"id1" : 1, "id2" : 2}

and the parameters for pointsalonglines algorithm, something like this :
param2 = {'INPUT':layer,
          'DISTANCE': QgsProperty.fromExpression(mydic["idfield"]),
          'START_OFFSET':0,
          'END_OFFSET':0,
          'OUTPUT':'test'}

where idfield is the name of an existing field in my layer, with values = "id1", "id2", ...
I tried many ways to write my expression but without any success.


texts: Order of KML Layers in Google Maps API V3

texts:
I'm working on a website that when finished will allow the user to turn on and off a number of different KML layers using checkboxes, the layers are parsed through geoXML3 and I want to know how to control the order which the layers appear.

Currently if I have 2 polygon layers (layer A and B) and I click the checkbox to load layer A then I click the checkbox to load layer B layer B will be on top as it should be. However if I then uncheck both boxes, and select B first and then A, B is still drawn on top.

I assume this is because I parse the data only the first time the boxes are checked and afterwards I just show or hide them, but is there a way to manage which layer is on top?

Below is an example of the function I use to parse the layer and then turn it on or off.

function boundaryToggle_FOS() {
if (boundaryIsLoaded_FOS == 0)
{
boundaryIsLoaded_FOS = 1;
    if (firstLoadFOS == 0)
        {
        kml_FOS.parse('FOS_Boundary_v03.kml');
        firstLoadFOS = 1;
        }
    else
        {
        kml_FOS.showDocument();
        }
}
else
{
    kml_FOS.hideDocument();
    boundaryIsLoaded_FOS = 0;
}}


texts: Loading raster from ArcSDE with MSSQL does not work

texts:
I have a problem loading rasters from ArcSDE using MSSQL. The rasters are displayed as polygons with the geometry "footprint" and in QGIS only the outer border of the grid is generated as a polygon. All other data are displayed correctly.

the import in QGIS looks like this:

And this is the layerinformation:


texts: Why is every SELECT via a geometry function null?

texts:
I'm using Spatialite (2.3.1, binary from the website) on Windows XP via Java ("Xerial Driver"). SQLite loads fine, also Spatialite seems to be loading successfully. I can even do SELECT-statements. But I can't do SELECT's on geometries. E.g.

SELECT link_id, SRID(GEOMETRY), AsText(GEOMETRY) FROM streets1

returns only null columns except for link_id. How can that be? By the way, I imported this data originally from MapInfo files using ogr2ogr:

ogr2ogr -clipsrc 8.6688 50.1018 8.6794 50.1082 -f "SQLite" streets1.SQLite ~/Desktop/G2AM10101EG2000MAADT/Streets1.TAB

Also I tried:

ogr2ogr -clipsrc 8.6688 50.1018 8.6794 50.1082 -f "SQLite" streets1.SQLite Streets1.TAB nlt GEOMETRY

This also gives the same result...

Thanks, Philip


texts:
I'm using tilestrata and mapnik server to host geotiff map.
I'm getting very strange behaviour:

Basically, rendered tile shifts across map in strange fashion.
The transformation from lat/long to XYZ is known and I implemented it correctly (works on other map).

One lead I can think of is the fact, that the map uses different units for pixel width and height (other maps I was using had square pixels).

EDIT: Converted map to have "square" pixels using gdalwarp and gdaltranslate -  bug is still on!

Any ideas what can cause such a strange behaviour?

Can I somehow set different aspect ratio of pixel in proj4 srs info?

Here is the gdalinfo from the file:

Driver: GTiff/GeoTIFF
Files: GIG_halda_jadwiga_orto_2cm.tif
Size is 40549, 35453
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (19.165103727574930,50.334287477469992)
Pixel Size = (0.000000280890000,-0.000000179799000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  19.1651037,  50.3342875) ( 19d 9'54.37"E, 50d20' 3.43"N)
Lower Left  (  19.1651037,  50.3279131) ( 19d 9'54.37"E, 50d19'40.49"N)
Upper Right (  19.1764935,  50.3342875) ( 19d10'35.38"E, 50d20' 3.43"N)
Lower Right (  19.1764935,  50.3279131) ( 19d10'35.38"E, 50d19'40.49"N)
Center      (  19.1707986,  50.3311003) ( 19d10'14.88"E, 50d19'51.96"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Overviews: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Overviews: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Overviews: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
  Mask Flags: PER_DATASET ALPHA
  Overviews of mask band: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha
  Overviews: 20275x17727, 10138x8864, 5069x4432, 2535x2216, 1268x1108, 634x554, 317x277, 159x139


texts: How to exclude parts of route based on street features?

texts:
I have a task. I have a route (coming from whatever, I am using Google Maps for now, but I can switch to OSM), and I want to exclude some parts of this route based on street features (like I don't need highways). I need to do it automatically on node.js server. Can you suggest any tools I can use on server side to solve this problem?
EDIT: I added some info, that I need to do it not manually, but automatically on node.js-based server.


texts: ArcGIS Javascript API: Force the fit option when zooming via rubberband selection with tiled map services

texts:
We're using ESRI's ArcGIS Javascript API with a tiled basemap and when the user zooms in using a rubberband selection via the navigation toolbar the map sometimes zooms in too far.  This is apparently expected behaviour (not sure why!) when using the API with a tiled map service.  However, we need to change this behaviour since this is a usability issue.

When calling map.setExtent() directly it's possible to change the default behaviour by specifying the fit option.  Eg. map.setExtent(extent, true).  Unfortunately, we can't find a way to specify this option when zooming using the Navigation toolbar.  The only way to fix this that we could find was by overriding setExtent() in a custom sub-class of esri.Map.  We ended up having to override the setLevel() method as well because just blindly setting the fit option to true inside our setExtent() method caused the setLevel() method to fail as the zoom level was sometimes overridden.

Functional workaround:

dojo.declare("MyMap", esri.Map, {
  levelling : false,

  setExtent : function(extent, fit) {
    var shouldfit = fit;
    if (fit === undefined && !this.levelling) {
      shouldfit = true;
    }
    this.inherited(arguments, [extent, shouldfit]);
  },

  setLevel : function(level) {
    this.levelling = true;
    try {
      this.inherited(arguments);
    }
    finally {
      this.levelling = false;
    }
  }
});

I'm not happy with the workaround above and feel we must be missing an obvious, simple way to change the zoom behaviour to use the fit option without resorting to the hack above.  Does anyone know of a better way?


texts: Grey Lines when batch processing a self build arcmap tool with model builder

texts:
I've used ArcMap 10.3 ModelBuilder for a simple model that iterates the raster data in a folder and writes the result of raster calculation to another specified folder.
Instead of the specified input folder name and output filename I rename them to Input and Output. Processing this within ModelBuilder works fine.
But as I have a few folders to process, I tried the batch processing mode which ended with grey lines for the output file names.

How can I modify my model in order to be able to process this simple model in batch mode?

Edit 1:
I'm trying to iterate the rasters in the folders (examp: E:/1/path1raw1/; E:/1/path2raw2/; E:/1/path3raw3/) and try to save them to another folder (examp E:/2/path1raw1/; E:/2/path2raw2/; E:/2/path3raw3/) after calculation by specifying the output as E:/2/path1raw1/%name%.tif. I think the problem is that it can not iterate the output folders based on input folders. But I can not find a solution yet.


texts: Data Interoperability FME - Python Exception : u'hasZ'

texts:
I have an issue with data interoperability using a feature service as a reader and a writer.

I'm reading from a feature service a table and using it to enrich it with data from an Access database and send back the updated version.

Two readers:

-ArcGIS Portal - Feature Service "A"

-Access Database

One Writer

-ArcGIS Portal - Feature Service "A"

I receive an error as follows:

  Consider using coordinate system 'EPSG:31700' instead of '_EPSG:31700_0'

  ArcGIS Portal Feature Service Writer: Getting item with ID '3a50c62127f8473fa7feb8be1be0c483'

  ArcGIS Portal Feature Service Writer: Requesting metadata for Feature Service at 'http://--------------.ro/server/rest/services/Hosted/Cornu_test/FeatureServer'

  ArcGIS Portal Feature Service Writer: Getting item with ID '3a50c62127f8473fa7feb8be1be0c483'

  ArcGIS Portal Feature Service Writer: Requesting metadata for Feature Service at 'http://--------------.ro/server/rest/services/Hosted/Cornu_test/FeatureServer'
  ArcGIS Portal Feature Service Writer: Found layer/table 'Person' with ID 3 for feature type 'Person'

  ArcGIS Portal Feature Service Writer: Requesting metadata for layer/table 3 of Feature Service at 'http://--------------.ro/server/rest/services/Hosted/Cornu_test/FeatureServer'

  Python Exception : u'hasZ'

  A fatal error has occurred. Check the logfile above for details

  A fatal error has occurred. Check the logfile above for details