texts: Buffer point dataset based on compared value

texts:
I am creating polygons from the dissolved buffers of points and need to control the dissolve process. In the model the first stage is to assign high/med/blackspot using the following

f2='ERROR'
def re_score(my_field,f2):

 if (my_field*1) < 14:
  f2 = 'Black'
 elif (my_field*1) < 24:
  f2 = 'Medium'
 elif (my_field*1) >= 24:
  f2 = "High"

Expression = re_score(!Rx!,f2)

and then use the following settings in buffer...

arcpy.Buffer_analysis(Signal_Strength_Added, Buffered_tower_range, Distance__value_or_field_, "FULL", "ROUND", "LIST", "Signal_Strength")

If distance between points is 30m then I buffer at 35m

I want to deal with overlap by defaulting to the worst case scenario, so if when dissolving the buffers are "med" and "blackspot", then it should be incorporated into the blackspot zone - I don't know how to specify this!

Any ideas?

The base points dataset is

Latitude    Longitude   Rx_dB  Best_unit
-17.06694    145.23158  -050.5  2
-17.06695    145.23297  -044.1  2

my_field in this case is the Rx_dB

The image below illustrates the issue...

_pts is the output point dataset with signal strength
_Lamb is the 30m buffer with dissolve by signal strength
_pts_30m_Buffer - is just the buffer without the dissolve

Areas around 1a should be red, 2a/2b area should be yellow and 3 is correct.

Each dataset is around 500,000 points and I am still trying to do a spatial join on a new dataset without the dissolve.


texts: Generate point layer (grid values) with Z coordinate of (TIN) interpolated raster in QGIS 3.16

texts:
1. Question With Additional Information:
I would like to generate a sparsely packed layer containing points labelled with the Z-column of an existing interpolated TIN raster
A sample photograph what is needed is attached below,

2. Background Information:
The following explains my process flow,
Step 1. Raw data in the form of .csv files; one of which contains
(i) depth measurements at various locations from, say, a water body, i.e., 3 columns of latitude, longitude, and altitude (altitude represents the Z-column described in the question), and the other.
(ii) boundary of the area covered previous explained in (i); again in latitude, longitude and Z-column form.
Step 2. Generate .shp files using OSGeo and GeoPandas - Python of items in Step 1 (along with the required CRS).
Step 3. Using TIN interpolation to generate raster of item described as (i) in Step 1.
Step 4. Use Clip Raster by Mask Layer to remove extents beyond the boundary described in (ii) in Step 1. The function, minimum bounding geometry is used to create the required mask layer from the boundary layer csv file.
3. Step of Interest Wherein the Question Asked Arises:
Step 5. I, now, use Pixels to Points to generate the point layer and further, add Labels->Single Labels-> Select the column containing the Z-column. However, the result obtained is very densely packed that the underlying layer is no longer visible. This is described in the image below,

How can I modify the density/ resolution/ the number of points to achieve a result close to that described in the Image in 1?


texts: How to batch "Layer save as"-process in QGIS?

texts: I want to save many (about 100) QGIS-layers at once. Is there an automatic way to batch the "save as"-process?


texts:
I have a csv-file and want to load it in the eVis plugin. The csv-file has a column cid and is sorted ascending.

sorting in csv-file:

I have pictures of a road and this pictures should be shown in the correct sorting.
How can I mange this in the plugin?

I tried different ways to import the csv file to QGIS, the plugin spredsheet layers, drag & drop and the delimeted textlayers, but got always the same result.

The first object in eVis but with the cid 1009


texts: Selecting one polygon for one point using ArcGIS Pro

texts:
I have a points shapefile and a polygon shapefile. I work in ArcGIS Pro.
I would like to select all the polygons which have one or more points in them, however, one condition should be considered, that if one point lies in two or more polygons which overlays, only one polygon is selected.
Please see picture for example. There are 3 polygons highlighted with orange circle with only one point in them. As the polygons overlay, the result I obtained from ArcGIS Summarize Within tool showed that each of these polygons consist of one point, which is wrong, cause there are not 3 points, but 1. There should be only one polygon with one point, the other one stays empty.

The result should look something like this. Orange points are selected polygons with one or more points. And pay attention to black line - these are polygons which overlays and shares a point, but only one of them are selected. Of course this could be done manually, but it's too time consuming and it gets hard when there are to many overlaying polygons, so I try to find a fast way to do this.

Here is a very similar question Selecting polygons which contain at least one point with spatial index about at least one and more points in polygon, but the solution selects all polygons with intersection, even if they overlay and share a common point.
I couldn't find the solution with "Join attributes by location" or "Select attributes by location". I'm probably missing something and this should be done easily.
Is there some fast and effective way to do this?


texts:
I'm triyng to use the python bindings of gdal to rasterize a raster using a shapefile.
The code I'm using is rather simple:

gdal.RasterizeLayer(target_raster, [1], input_layer ,burn_values=[-9999],options=['INVERSE=TRUE'])

The code runs and doesn't produce any errors, but it is completely ignoring the "inverse" option meaning that it is removing the part of the raster I want to keep.


texts: Weighted polygon-feature-classes combination/overlay analysis using ArcGIS for Desktop?

texts:
I have several (sometimes, not always overlapping) Polygon Feature Classes which are all weighted/ranked through a field from 1 to 10 (10 is best). I would like to combine/overlay these layers them, but also use their rank information. Therefore not only the overlapping/intersected features should be shown in the result. A classification of the "best" places with the highest combined ranking should be shown in the result.

Is there a way to do that using vector tools?


texts: Setting default value for date columns with FeatureClassToNumPyArray in ArcGIS Pro?

texts:
I am reading a feature class into a pandas dataframe using:

df = DataFrame(
    arcpy.da.FeatureClassToNumPyArray(
        in_table=feature_class,
        field_names=field_list,
        skip_nulls=False,
        null_value=np.NaN
    )
)

Where feature_class is the name of the feature class and field_list is the list of fields to include.

Setting skip_nulls=False ensures the floats are converted to np.NAN and setting null_value=np.NAN does the same for integers. However, one of the columns has a date type, which causes the following error:

  ValueError: year 0 is out of range

Is there any way to set a specific default for the date type?


texts: Strip Map Index Tool: output has no geometry

texts:
I'm having a problem creating a grid index in ArcMap. I'm encountering the same problem with both the Strip Map Index Features (first choice) and Grid Index tools.

Steps:

Open the Strip Map Index Features tool in the Cartography toolbox.
Specify a line feature as input and set output to default geodatabase.

Result:

--Output features were generated, but have no geometry. "Shape_Length" and "Shape_Area" (as well as "TopPage" and "BottomPage") fields in attribute table contain null values. Symbology is set, but there's nothing to display.

I get a comparable result when I attempt to use the Grid Index Features tool.

System:
ArcGIS 10.2.1 (Basic license) | Windows 10 (64-bit)

I was able to work around this problem by creating a whole new MXD and adding all my data back in by hand. Not a definitive solution by any means.


texts:
Working in consulting, I found that sending clients PDF file with my embedded maps is one of the best ways to achieve a wow-effect from them.

The latest pdf file formats all support the layers feature (which arcmap makes use of); however they also are able to support javascript.

Does anyone have ideas to use the power of javascript (or anything else) to make my maps within my pdf files more dynamic?

Which tools do you recommend? Any extensions to esri products?


texts:
I have created a PostGIS DB using osm2po and I am using pgrouting to calculate the routes.

The application I am creating should find routes for cyclists so I want the route finder to prefer roads with cycle lanes present and to never plot a route containing motorways.

Looking through the config files I see: wrt.finalMask = car along with other tags such wtr.tag.highway.motorway = 1, 11, 120, car

I was wondering are these the parameters that I should tweak to accomplish this?

Would I just remove tags such as motorway etc.?

Thanks


texts:
From Mapbox's documentation, I should get a GeoJson feature object when I call queryRenderedFeatures method.

I run it from a click event without any filter:

this.map.on("click",
    function(e) {
        var clickedCluster
        this.map.queryRenderedFeatures(e.point).forEach(function(element) {
            console.log(elem)
        })
    })

On my console (see screenshot below), I got feature's properties of the point I click on the map, I got layer's information, but I don't get its geometry. How can I get it?

As far as I understand Mapbox's example, I should obtain it
https://www.mapbox.com/mapbox-gl-js/example/queryrenderedfeatures/


texts: How to change projection to lat long from eastings and northings, QGIS

texts: I'm using QGIS and I'm trying to get the coordinate system into lat/long instead of eastings/northings. I think I'm choosing the correct coordinate system, but my "Extent" is listed in eastings/northings. I'm new to QGIS so I don't know what I should do next. Is there a re-project tool that I'm missing? I've attached an image for reference:


texts: OSM Buildings - problem adding .obj file

texts:
I've started playing with OSM Buildings, have modified an example.
I am unable to display a custom object I'm downloading using map.addOBJ
The script downloads the .obj and .mtl files from the webserver.
In the console, I see a "Uncaught TypeError: Cannot read property 'color' of undefined
at processOBJ (blob:null/afde6e25-7336-40a4-8805-46fbdeb82782:1)
at XMLHttpRequest.r.onreadystatechange"
Here's the code:
<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="width=device-width, initial-scale=1.0, maximum-scale=1.0, user-scalable=no">

  <link href="https://cdn.osmbuildings.org/4.0.0/OSMBuildings.css" rel="stylesheet">

What to do?


texts: Is Terrain and Topographic Correction the same for Landsat8?

texts:
Landsat 8 L1 products are "Terrain Corrected",

that said it means that are they already "topographically corrected", like both words are synonyms or you need to apply correction techniques as well? (eg i.topo.corr for Grass)


texts: Unable to Edit Shapefiles - New Installation

texts: I recently installed 10.6 and everything seems to be working great except I can not edit. It began every time I tried to edit a shapefile or feature class- the map would crash even with no basemaps or other data, so I added the editor toolbar and it seems to be "missing" something which explains why I am having this issue.


texts: Converting raster to polygons in ModelBuilder using FOR iterator gives ERROR 000622

texts:
I am having trouble with a section of a larger ModelBuilder routine. I have extracted all but the offending tools to see if I could get those to work. Basically, it takes a raster DEM and converts it to an integer raster based upon a cutoff value. It then converts the integer raster to polygons using the Raster to Polygon tool. It is at this point that it fails.
The model employs a FOR iterator to go through the DEM at increasing elevations of 10. If I run it through just once (beginning and ending elevations are equal), no problem. But if I iterate it two or more times, I get this error:

ERROR 000622: Failed to execute (RasterToPolygon). Parameters are not
valid.
ERROR 000861: Field:  is not valid.
Failed to execute (Raster to Polygon).

A couple of things to add:

If I execute it from the Edit menu using Run, it works fine and I get multiple outputs.
There are no spaces in the paths to the files. Changed this to the default .gdb with the same problem.
The output polygons do not approach the 2.1 gb limit for features I noted in the ESRI documentation.
I am running it in foreground.
There is definitely data at each stage of the iteration. I tested this by running each level as a single pass.
I am using %value% to change the intermediate file names for each pass.

I am at a loss as to what I am doing wrong.


texts: Clip SpatialPolygonsDataFrame and preserve data in R

texts:
I'm trying to clip a SpatialPolygonsDataFrame and keep it as a SpatialPolygonsDataFrame. However, when I use gIntersect, it converts it to a SpatialPolygons and I lose all the data. I'm pretty new at working with spatial data in R, so I'm probably missing something simple.

As an exercise, I'm trying to cut off the outlying islands of Hawaii.

library(tigris)
st<-tigris::states(cb=T)
hawaii<-subset(st,STATEFP=="15")
class(hawaii)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

hawaii2<-gIntersection(as(extent(-178.3347+18, -154.8068, 18.91036, 28.40212), "SpatialPolygons"), hawaii, byid = TRUE, drop_lower_td = T)
class(hawaii2)
[1] "SpatialPolygons"
attr(,"package")
[1] "sp"


texts: Converting Web Mercator (3857) to NAD 1983 StatePlane Texas Central FIPS 4203 (2277) in ArcMap?

texts:
I am trying to convert Web Mercator (3857) to NAD 1983 StatePlane Texas Central FIPS 4203 (2277) in ArcMap but whenever I use the Project tool (Data management) it says there are no transformations available.

What is the correct way to do this transformation using ArcMap 10.5?


texts:
I am trying to write two layers into geopackage without success.

from qgis.PyQt.QtCore import *

def drawLay ( namee, x1, y1 ) :
    lay = QgsVectorLayer('Point?crs=epsg:2039', namee , 'memory')
    pr = lay.dataProvider()
    lay.startEditing()
    fet = QgsFeature()
    fet.setGeometry( QgsGeometry.fromPointXY(QgsPointXY(x1,y1)) )
    pr.addFeatures( [ fet ] )
    lay.commitChanges()
    lay.updateExtents()
    QgsProject.instance().addMapLayer(lay)

# create two layers
drawLay ( 'aa', 150000, 400000)
drawLay ( 'bb', 150100, 400000)

#write this layers into geoPackage
GpName = r'c:\temp\gpk_file.gpkg'
options = QgsVectorFileWriter.SaveVectorOptions()
options.actionOnExistingFile = QgsVectorFileWriter.CreateOrOverwriteLayer
options.EditionCapability = QgsVectorFileWriter.CanAddNewLayer
for layer in  iface.mapCanvas().layers() :
    namee =layer.name()
    print (namee)
    options.layerName = namee
    _writer = QgsVectorFileWriter.writeAsVectorFormat(layer, GpName, options )


texts: How to Split a Concave Polygon Into Convex Polygons in Arcpy?

texts:
I'm looking for a tool or algorithm to detect concave polygons and split them into convex polygons. Like explained in the picture, the blue polygon is split into A and B polygons

I'm using Arcpy under Arcgis 10.1


texts: QGIS Installer switch to set target directory

texts: Can anyone please describe the command line switches that can be used when installing using the stand alone installer for ms windows. As far as I know one could use /S for a silent but I'm pretty sure there are other options as well such as choice of target directory but I just can't find them.


texts: Enhance the color associated with a band od an aerial photo?

texts:
I am using QGIS 2 on Windows 7 and have 4-band aerial photos loaded in to a project. I have switched the bands in order to make the photos display as color infrared. Does anyone know if you can enhance the color of the 4th band to make the red show better? I have a slight color blindness issue that hinders me seeing some of the red.

FYI... I am attempting to do a desktop analysis for a habitat assessment.


texts:
I have two digital elevation model rasters. I calculated the difference of these two rasters using raster calculator and got a single output raster. Now I want to highlight areas where the surface has dropped or eroded.

Is there a way in ArcPro to select all of the pixels in the new raster that have negative values?


texts: pyqgis processing runalg release input in windows

texts:
I am having problems with deleting files on disk after using "processing.runalg()" in a plugin I created when running it in windows. I have no problem in linux.

I found an issue regarding the same topic at https://hub.qgis.org/issues/12502, and was wondering if the issue in had been fixed in qgis 2.18.

The code is simple as follows, where outName is the path to the output chosen by the user, and tarNames is a list of file names to be processed:

import processing
import sys, gdal_merge as gm

rlayer = QgsRasterLayer(os.path.join(out,inputName), "input")

args = ['-of', "GTiff", '-n', "0", '-a_nodata', "0", '-separate',
    '-ot', "UInt32",'-o', out+"/"+tmpName]
args.extend([out+"/"+t for t in tarNames])

sys.argv = args
gm.main()

processing.runalg("gdalogr:warpreproject", {"INPUT":os.path.join(out,tmpName),
            "SOURCE_SRS":"EPSG:"+str(rlayer.crs().postgisSrid()),
            "DEST_SRS":"EPSG:32724", "NO_DATA":0, "METHOD":0,
            "RTYPE":3, "OUTPUT":os.path.join(out,outName)})

os.remove(os.path.join(out,tmpName))

os.remove fails, saying that the file is not available. As I mentioned, I don't have this issue in linux.

I tried to implement the workaround proposed in Releasing PyQGIS file locks?, but it does not solve the problem.

The only solution I have for the moment is to create a bunch of temp files with unique names and have the user delete them by hand... less than optimal.
The temp file is created on disk by gdal_merge.py, I did not find an option to create it in memory. If anybody has a solution for that, I am very interested. I understand this is a separate question, and I am happy to post it if there is a solution.


texts: Viewing Rasters on Phone/Tablet with GPS Location?

texts:
I have a Geotiff like the one below, and I'd like to be able to use it as a map I can view with my location on a table or smartphone. This can be via Google Maps or any other app.
That is, I can use this app and my phone's GPS to see where I am on this map. That way, I can send/share the maps with clients and they can view it on their property.
I don't need to be able to capture or work with the data in anyway. I just want an overlay with GPS location usable by a smartphone or tablet.
I have tried this with Google Maps, but as it can only work with vectors and I can't seem to get anything close to this image quality with a Vector it doesn't seem to work. And because I want clients to be able to use this, a $75 app like Guru Maps (that every client would have to buy) isn't a solution either.
I feel like there is an easy answer I haven't found yet.


texts: python.multiprocessing and “FATAL ERROR (INFADI) MISSING DIRECTORY”

texts:
I have a question which is closely related to this old question here: python.multiprocessing and "FATAL ERROR (INFADI) MISSING DIRECTORY"

I have also looked at this, which seems to relate to something similar: Multiprocessing Errors - ArcGIS implementation

just like in the linked question, I run into the error: "FATAL ERROR (INFADI) MISSING DIRECTORY" at random times, but never the BND-error

I have multiple folders with thousands of shapefiles in each, where I want to convert two of the shapefile fields to two different rasterlayers. the shapefiles were converted from csvs via multiprocessing, with no issues whatsoever.

I have tried most of the workarounds from the previous question, such as creating temporary and unique scratchworkspaces and workspaces, check that the infile exists and create a string in memory and flush it, but to no avail. I have also tried splitting the workload, so that it will first write all rasters for one field, then do all the rasters for the second field.

so my question is, would anybody know of another potential workaround? I was thinking that it might be possible to somehow make sure that when a new raster is written to the folder, that no other processes try to access the space, but I'm unsure how to code it.

relevant code excerpts below:

import arcpy
import os
import sys
import arcgisscripting
import multiprocessing
from functools import partial

def multiproc(infile,path,savepath1,savepath2):

 try:
        arcpy.Exists(infile)

        # Set local variables
        inLayer = path + "sample.xml"
        fieldList = arcpy.ListFields(infile[:-4]+".dbf", "field1")
        inData = infile + " field1"
        outLayer = infile[:-4]

        # Check out the ArcGIS Geostatistical Analyst extension license
        arcpy.CheckOutExtension("GeoStats")
        # Execute CreateGeostatisticalLayer
        arcpy.GACreateGeostatisticalLayer_ga(inLayer, inData, outLayer)
        arcpy.GALayerToGrid_ga(outLayer, savepath1 + "/" + infile[:-4])
        arcpy.Delete_management(outLayer)

        #starting second raster creation process
        arcpy.Exists(infile)
        fieldList = arcpy.ListFields(infile[:-4]+".dbf", "field2")
        inData = infile + " field2"
        outLayer = infile[:-4]

        # Check out the ArcGIS Geostatistical Analyst extension license
        arcpy.CheckOutExtension("GeoStats")
        # Execute CreateGeostatisticalLayer
        arcpy.GACreateGeostatisticalLayer_ga(inLayer, inData, outLayer)
        arcpy.GALayerToGrid_ga(outLayer, savepath2 + "/" + infile[:-4])
        arcpy.Delete_management(outLayer)

    except:
        import logging
        LOG_FILENAME = 'c:/error.txt'
        logging.basicConfig(filename=LOG_FILENAME,level=logging.DEBUG,)

        logging.debug(infile)

        logging.exception('Got exception on main handler')
        pass

if __name__=='__main__':
# Set environment settings
#setting paths for the directories

##Finding paths and list of files

#multiprocessing part
    Pool=multiprocessing.Pool
    workers=multiprocessing.cpu_count()
    pool=Pool(processes=workers)

    partial_multiproc=partial(multiproc,path=path,savepath1=savepath1,savepath2=savepath2)
    pool.map(partial_multiproc,filelist)
    pool.close()
    pool.join()

after having moved the data and tried on an internal harddrive (the previous attempts have been on an external HD), I also got the BND-error on rare occassions, whereas the INFADI is much more prevalent.

I have been unable to find a full solution so far, but a potential workaround lies in creating temp-folders that are deleted with arcpy at the end of the multiproc-function. even though the error persists, the faulty folders are - as far as I can tell - not deleted, meaning that at least it is possible to know which rasters were not completed and redo these later.


texts: trouble installing the qgis2kml plugin (or, how to install plugins not available in the Plugin Manager from source)

texts:
I'm trying to install the qgis2kml plugin for qgis. Its not showing up in the Plugin Manager search and the only qgis "repositories" I can find are websites for downloading source code.

Which is fine, except for that I don't know how to install python source into qgis, and there doesn't appear to be a readme available in the plugin source code.

I could solve my problem, installing the qgis2kml plugin, in three different ways, based off what I've described so far:
- by installing directly from the source code using some low level apis or some such thing
- by finding a repository aside from the Official XML repository that has qgis2kml

OK so that's two. I'm wanting to export a heatmap from qgis that is colorizing over 50 different shapes in 9 different shades of a segmented continum from green to red. This would take so long to do I simply wouldn't do it and it would be amazing to export heatmaps as translucent kml overlays as I know can be done in kml.


texts: Why might Spatial Join appear to produce incorrect values?

texts:
I'm trying to do a spatial join of a point file (created from a raster) to a polygon (target), where I am joining the GRID_CODE value from the point to the polygon. I've tried a variety of permutations, one-to-many, many-to-one, changing the order of the input files, putting in a mode and mean condition for field_map, checking the field properties (both Double), trying a multi-part to single-part conversion (data management) for both inputs, but the resulting joined field is always populated with incorrect data (as verified by looking at the point values).

For example, here is a screenshot of the correct values from the point file, highlighting the 2 (green)

Here is a screenshot of the spatial join output highlighting that same 2 value, which I would expect to be in the same location.

What is additionally bizarre, and might be a clue as to what is going on, is that the same field (GRID_CODE) has a different number of values (4 in the original, and 5 in the output

Any suggestions would be much appreciated!


texts: Figuring out whether a feature is visible or not in pyqgis

texts:
I have a vector layer (lines) with a numeric (integer) attribute, let's call it Status. This layer is styled according to the attribute using a categorized style:

So far, so good. Now I'm trying to figure out a way to get the visibility status (i.e., I click the checkbox on the left) of a feature in this layer using pyqgis. Here's what I have so far:

It's possible to get the individual QgsRendererCategoryV2 instances by using layer.rendererV2().categories(). This categories, however, don't seem to contain the visibility information, but only the values and colors (or I'm missing something).
The individiual symbols (e.g. with layer.rendererV2().symbols()) don't help me either.

Ideally I would like a solution like feature.isCurrentlyVisible(), but I'm prepared to having to look up the attribute in the table myself.


texts: Converting Text Field to Numeric in ArcGIS Pro?

texts:
I'm trying to convert text to numeric, but am unable to in ArcGIS Pro.

Is there a Python function?

I tried int(!FieldName!) but I think NULL values are messing it up?


texts: Check if lat/long point is within a set of polygons using google maps

texts: I know how to do the above using Esri (Query task on shapefile), but can this also be done using Google Maps?  Do I have to query each polygon or is there one method to query a set of polygons?


texts: How to automate creation of PNG maps using Python?

texts: I have a project with 50 vector layers and I have done custom styling for them. Now I want to make around 30 png files with different xmin, ymin, xmax, ymax. The output size (png) should be width 139.7 and height 152.4. Every quarter I have to run this exercise manually around 3-4 times. Is it possible to automate this?


texts: To consume geojson data from a REST HAL CORS webservice, in client side i'm using ExtJs 5, but to map data having geometries i'm stuck between choosing geoExt2 which has a branch based on extjs5 and openlayers 2.3, or working directly with extjs5 and openlayers 3, for consuming data till the moment i'm using rest proxy in extjs, but geoext proxy seems not to be agile because it's based on OpenLayers.Protocol.HTTP which do not support REST, is there any suggestions how to use extjs5 and ol3 with rest.


texts: Leaflet: resetStyle function doesn't work correctly

texts:
I made a map adapted from the leaflet tutorial for an interactive choropleth map. Highlight function is working correctly, but nothing happen when the mouse is out of a polygon: resetStyle doesn't work on my map!

This is my code:

var val = new L.layerGroup ()
val.addTo(map);

function valleescoul (feature) {
return {color: 'purple',
        fillColor: 'blue',
        weight: 2,
        opacity: 1,
        fillOpacity: 0.3};
   }

function highlightFeature(e) {
    var layer = e.target;

    layer.setStyle({
        fillOpacity: 0.7
    });

    if (!L.Browser.ie && !L.Browser.opera && !L.Browser.edge) {
        layer.bringToFront();
    }
}

function resetHighlight(e) {
    layer.resetStyle(e.target);
}

function zoomToFeature(e) {
    map.fitBounds(e.target.getBounds());
}

function survol(feature, layer) {
    layer.bindTooltip(feature.properties.name, {className: "labelsmall", noHide: false, sticky: true, direction: 'center'});
    layer.on({
        mouseover: highlightFeature,
        mouseout: resetHighlight,
        click: zoomToFeature
    });
}

  L.geoJson(vallees, {style: valleescoul, onEachFeature: survol}).addTo(val);

Where did I make a mistake?


texts: Fastest way of adding dictionary values to attribute table with PyQGIS3

texts:
I have a Python dictionary with around 450,000 keys and respective values. I want to join this dictionary to an attribute table of a vector dataset with roughly the same amount of features.

First I created a new field where the values are to be stored and I saved its index in a variable. I also saved the index of the field through which I want to join the values of the keys.

layer_prov = layer.dataProvider()
layer_prov.addAttributes([QgsField("new_field", QVariant.Double)])
layer.updateFields()
new_field_idx = layer.fields().indexFromName('new_field')
join_field_idx = layer.fields().indexFromName('join_field')

layeris the vector dataset that I want to join the information from the dictionary to. dict is the dictionary.

for key, value in dict.items():
layer.startEditing()
for feature in layer.getFeatures():
    attrs = feature.attributes()
    join_value = (attrs[join_field_idx])
    join_value = str(join_value)
    if join_value == key:
        feature[new_field_idx]=value
layer.commitChanges()

Is there any way of doing this faster? For every entry in the dictionary this model needs to go through all the features of the vector dataset (450,000 x 450,000 = 202,500,000,000 comparisons).


texts: Creating ArcGIS Raster Catalog with date/time fields from directory of NITF images?

texts:
I have a directory of NITF files, each with the idatim (timestamp) field set.  In the ArcMap properties, it shows up as "acquisition date" in the Key-Metadata table.  I would like to be able to create a raster catalog, mosaic, or something which will enable me to enable time support in ArcMap so I can create a time-lapse video.  The date field should be using this acquisition date.
I am getting familiar with ArcPy if that is a recommended solution. I have read through the ArcGIS documentation on using the time slider, how time is represented in datasets, and how to enable time on different products.  The key issue I am detecting is that when loading NITF rasters into a geodatabase, they do not set the key-metadata field as the date_time attribute.  I imagine there is a clever way of doing this, however I have not found it yet.  If this is not possible, what part of the NITF raster should provide the time.  I can update the NITF metadata myself if needed.
Some research
http://resources.arcgis.com/en/help/main/10.1/index.html#//005z00000004000000
http://resources.arcgis.com/en/help/main/10.1/index.html#//005z0000000p000000
http://resources.arcgis.com/en/help/main/10.1/index.html#//005z0000000q000000


texts: Area Calculation: Different result between GEE and ArcGIS

texts:
Can anyone explain why I got different area calculation result between GEE and ArcGIS software?
This my GEE code:

for(var x = 1; x < 5; x++) {

  var cal = image_class.eq(x).multiply(ee.Image.pixelArea())
                             .reduceRegion(ee.Reducer.sum(),
                              aoi, 20, null, null, false, 1e13, 1)

  print('ID' + x ,cal,ee.Number(cal.values().get(0)).divide(1e6).multiply(100))

}

These the results (hectares):

GEE:

ID: 1 = 13,637
ID: 2 = 30,314
ID: 3 = 5,663
ID: 4 = 648

GIS SOFWARE (PROJECTION: WGS_1984_UTM_Zone_49S):

ID: 1 = 15,923
ID: 2 = 27,678
ID: 3 = 6,159
ID: 4 = 579


texts:
I have this query where I'm trying to two join tables with st_intersects. There will be different geometry types in the workarea column but the point_geom will always be a point.

with cte1 as (
          select t.number,
             upper(t.workstreetname) as workstreetname,
                     upper(t.workcounty) as workcounty,
                     upper(t.workplace) as workplace,
                     upper(t.workarea) as workarea
          from schema.ticket t
              )
cte2 as (select distinct number, workstreetname, workcounty, workplace, workarea from cte1)
select c.number,
       c.workarea as workarea,
       a.geom as point_geom,
from cte2 c
join schema.addresses a on st_intersects(a.geom,c.workarea);

But I am confusing it somewhere because I get this error:

ERROR:  function st_intersects(geometry, text) is not unique
LINE 27: join schema.addresses a on st_intersects(a.geom,c.wo...
                                    ^
HINT:  Could not choose a best candidate function. You might need to add explicit type casts.
********** Error **********

ERROR: function st_intersects(geometry, text) is not unique
SQL state: 42725
Hint: Could not choose a best candidate function. You might need to add explicit type casts.
Character: 1246

Any suggestions on where I am going wrong on this?


texts: Using Python Add-In to delete/remove ArcMap layers that are turned off?

texts:
I have made a python addin button which removes layers that are turned off in an MXD.  The strange thing is that I get different results depending on whether I run the python commands in the python window vs my actual addin button.

In my addin class, the code looks like this:

def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        layerList = arcpy.mapping.ListLayers(mxd)
        for layer in layerList:
            if layer.visible == False:
                arcpy.Delete_management(layer)
        arcpy.RefreshTOC()

So, this works fine in the python window but doesn't do anything in the button.  I assume because it's trying to delete a layer object instead of the name.  I have made a workaround which changes the following line:

arcpy.Delete_management(layer.name)

Now this actually works.  But there is a caveat, if I have multiple layers with the same name but different visibility settings, it may delete the one that is visible.  My second workaround is to change the layer name before I delete it.

My final code looks like this:

def onClick(self):
        mxd = arcpy.mapping.MapDocument("CURRENT")
        layerList = arcpy.mapping.ListLayers(mxd)
        for layer in layerList:
            if layer.visible == False:
                layer.name += '_'
                arcpy.Delete_management(layer.name)
        arcpy.RefreshTOC()

By changing the name of the layer, I'm ensuring that I only delete invisible layers.  But I'm wondering if this is the best way to go about doing this.

My question is, is there a better way to remove layers that are not visible from an MXD?  Especially when they have the same name?


texts: Checking for values in one column and do something in another using ArcGIS Field Calculator?

texts:
I have two columns: "Street1" and "Verified".

I need to check all rows in Street1 to see if they are of string "0" or null.
I've attempted to come up with something in Field Calculator but it is returning an error saying the field [Verified] is not nullable


texts: Assigning legend items to specific column using Data Driven Pages

texts:
When using Data Driven pages, I'd like to be able to assign legend items to specific columns. I have each of the 3 top layers set to "only show classes that are visible in the current map extent." When all layers are present, the legend looks like this:

When the 2nd layer (SADC Easements) is not visible, it does not show the layer, but also doesn't push the "place layer in a new column" property to the next layer, so the legend shows up as:

This is a problem when trying to fit the legend in the space available for it and still have it a decent sized font for legibility.

I would like to have 2 columns regardless of if one of the legend items is not on the map for the current data driven page, like this:

Without having this when all the layers are on a data driven page:


texts: ArcMAP / R: Convert raster pixels to Long and Lat

texts:
I have downloaded a raster, which has the following information.

Item    Data
Date    2015-01-03 00:00:00
Product cold partial sst
Satellite   NOAA-18
Direction   Unknown
Dimension   1596 x 1617
BPP 8
Latitude range  47 — 62.999108
Longditude range    -15 — 13
Real → DN   dn = (rw - -3) / 0.1
DN → real   rw = 0.1 * dn + -3
DN range    0 — 255
Pixel dimensions    1.100149 x 1.100149

On Arcmap it shows the pixel X/Y co-ordinates. How can I convert this to Long and Lat?

From the FAQ of the site that I downloaded the raster from the following information is provided. link

How do I convert image coordinates to latitude/longitude?
For converting from x/y pixel co-ordinates on a Mercator projection image to latitude/longitude (lat/lon) positions use:

rows - number of rows in image
cols - number of columns in image
minlon/maxlon - minimum/maximum longitude of image (in decimal degrees)
minlat/maxlat - minimum/maximum latitude of image (in decimal degrees)
Note that these are the coordinates of the CENTRE of the corner pixels.
DEGTORAD - conversion from degrees to radians (PI/180.0)
RADTODEG - conversion from radians to degrees to radians
ln - natural log
Calculate longitude.

lonfract = x / (cols - 1)
lon = minlon + (lonfract * (maxlon - minlon))
Calculate latitude.

latfract = 1.0 - (y / (rows - 1))
Ymin = ln (tan (DEGTORAD * (45.0 + (minlat / 2.0))))
Ymax = ln (tan (DEGTORAD * (45.0 + (maxlat / 2.0))))
Yint = Ymin + (latfract * (Ymax - Ymin))
lat = 2.0 * (RADTODEG * (arctan (exp (Yint))) - 45.0)

The maximum/minimum latitude/longitude values and sizes are given when the images are viewed.

Would it possible to run the above code in ArcMAP or R (raster package)?


texts: ArcGIS Server - stopserver.sh vs. Admin Stop

texts:
I'm curious to the underlying "cleaning" differences between restarting the entire ArcGIS Server process with stopserver.sh vs. doing a stop through the arcgis/admin console.

http://server:port/arcgis/admin/machines/<machine>/stop

Obviously stopserver.sh shuts down the entire process, and admin-stop will keep up the admin/rest APIs and Manager.

I've noticed that the above command shuts down all services and arcsoc.exe processes, but does it keep open states in the background that get resumed on startup, or would this be the equivalent of running stopserver.sh in regards to application server maintenance?

The problem that I am running into is that over time that ArcGIS Server does a poor job of recycling resources, eventually causing issues with services. I'd like to run a scheduled reboot, but currently have a clustered configuration and this would require running stopserver.sh separately on each system.

I would rather create something that is system-independent (In case I migrate), and that works through the web adapter.


texts:
I have a hyperspectral image on which I have performed PCA and now intend on using the output PCA components as an input into a classification.

This is the first time I have done PCA, am i correct in thinking that I need to use the output of the Inverse PCA as input into the classification?

The Inverse PCA effectively recreates the original image from the PCA components and therefore has spectral properties, unlike the first outputs from the original PCA.

I presume people have done something similar - am I correct in using the Inverse PCA outputs for classification?


texts: ArcGIS ERROR 000824: The tool is not licensed in arcpy

texts:
I try to run the Zonal Statistic as a Table through arcpy and it gives me:
ERROR 000824: The tool is not licensed.

Here is the doc for the function: http://resources.arcgis.com/en/help/main/10.1/index.html#//009z000000w8000000

I have tried including the line where it checks for my spatial analyst extension:

arcpy.CheckOutExtension("Spatial")

but that just leaves the script running for a long time, way longer than it should (10 minutes+ on a 250Mb raster!)

I was hoping anyone else has a workaround/fix? Thank you for your time!

I have installed 10.1 Service Pack 1. I have checked and enabled all the extensions. Furthermore, I have the advanced ArcInfo version.


texts:
As we know esri shared their library which makes it possible to work with FileGDB format without using ArcObjects. So the question is what kind of applications could be created in order to help GIS specialist to solve their tasks related with this format? For instance, I have an idea to create an opensource tool(GPL) which converts FileGDB to Shapefile (or other format) and backwards. Will it be in the demand? or perhaps somebody have already developed such tool?
I will be grateful for any ideas about that (as well as collaboration ideas).


texts: Is there a way to programatically create a Mobile Project for ArcGIS for Mobile?

texts:
I need to build a tablet application using ArcGIS Runtime for Windows Mobile. Is there a way to create a mobile project usable by the ArcGIS for Mobile app without having to use the Mobile Project Center?

I need to be able, with a single click, to generate a shapefile (or pgdb), convert it to a mobile cache and be able to generate a mobile project out of it.


texts:
I am calculating WGS84 coordinates of points in an aerial image, given the WGS84 coordinates of the center of the image and the ground distance width and height (calculated from altitude and lens field of view).  Images cover a ground area of no more than 150x100m.  Thus far, I have neglected the effects of any coordinate system, and linearly calculated coordinates from the center coordinate.

Some inaccuracy in the coordinates has been noticed, particularly near the edges of images, and I am wondering if my neglect of coordinate system and projection is at fault.  If I were to take it into account, what "projection" would my images be in?  Could I simply convert to UTM Zone X for the computation, and convert the final result back to WGS84, or is UTM not an accurate representation of an image either?


texts: Is there an ArcGIS to PostGIS function cross-reference table?

texts:
I was looking back at the presentation called "PostGIS Tips for Power Users" by Paul Ramsey and got to thinking... Has anyone ever seen an ArcGIS Geoprocessing Tool to PostGIS function cross-reference table?  In particular, pages 21-27.  For ArcGIS users attempting to transition to PostGIS, I would think this could be very helpful.

I know of the following resources, but not aware of anything to help users learn the right terminology in PostGIS from what they already know in ArcGIS.

PostGIS 2.0 Cheatsheet
PostGIS Manual

For example, I was thinking something like this:


texts: How can I find credit consumption of my planet API key?

texts:
I'm using Planet API and I'll be making numerous requests every day. So the answers to this question will be a great help.

How to find credit consumption of API calls made on planet API?
Is there any section in Planet.com?
Will I be charged again if made the same API request again?

same API request meaning - Same input (item_id,item_type) to the same API(search imagery API, or download imagery API). Exact same request.


texts: How to add time in geoserver layer with gsconfig

texts:
Find below the code I use to add a raster to geoserver through python3 with gsconfig-py3 (1.0.7)

What I do

from geoserver import catalog
from geoserver.util import shapefile_and_friends
from geoserver.support import DimensionInfo

ws_name = "hiphen_iota"
store_name = "fda_2017_site03_ndvi_20170809"
style_tif_name = "hiphen_iota:ndvi_fda"

layer_path = 'U:/jeremy/FDA/Ortho_Indice/NDVI_Site03_20170809.tif'
layer_time = "2017-08-09T12:00:00.0Z"
###
# connect remotely to geoserver
cat = catalog.Catalog("http://##############:8080/geoserver/rest/")
# create / connect a new workspace
#ws = cat.create_workspace('python_test')
ws = cat.get_workspace(ws_name)
# create a new store with a shapefile
ft = cat.create_coveragestore(store_name, workspace=ws, data=layer_path )
# add a style
ft = cat.get_layer(store_name)
ft.default_style = cat.get_style(style_tif_name)
cat.save(ft)

Then I activate time dimension:

# add time dimension
timeInfo = DimensionInfo("time", "true", "LIST", None, "ISO8601", None)
ft.metadata = ({'time': timeInfo})

Afterwards, I found nowhere how to set time for my layer, I expect something like this:

What I want

ft = cat.get_layer(store_name)
ft.set_time("2017-08-09T12:00:00.0Z")
cat.save(ft)

My final purpose is to automatically add raster/layer to geoserver and use this: leaflet time dimension


texts: Controlling messages that are pushed to QgsMessageBar()

texts:
I've got a QGIS processing script which displays progress in the message bar, using the code in the PyQGIS cookbook.

It displays properly at the beginning of processing, but my code calls some other processing algorithms from Grass, which also push messages to the message bar and override my progress bar. Is it possible to either set the default duration of these messages to something small like a second, or block the message entirely?

I tried an approach that connects the widgetAdded slot to a function to re-add the progress bar, but it's causing QGIS to crash. Here is the code:

def putBarBack():
    iface.messageBar().widgetAdded.disconnect()
    iface.messageBar().pushWidget(progressMessageBar, iface.messageBar().INFO)
    iface.messageBar().widgetAdded.connect(putBarBack)

def run():
    #set up progress bar
    global progressMessageBar
    progressMessageBar = iface.messageBar().createMessage("Processing...")
    bar = QtGui.QProgressBar()
    bar.setRange(0,100)
    bar.setAlignment(Qt.AlignLeft|Qt.AlignVCenter)
    progressMessageBar.layout().addWidget(bar)
    iface.messageBar().pushWidget(progressMessageBar, iface.messageBar().INFO)

    iface.messageBar().widgetAdded.connect(putBarBack)

    # do some more processing here...

    bar.setValue(25)

    processing.runalg(grass7:r.watershed)

    bar.setValue(75)

    #do some more processing

    bar.setValue(100)


texts:
I have two raster layers: One with 56 m resolution and the other has 30 m resolution. Usually, I just use the 'Combine' command in ArcGIS to combine them, getting a single value for each unique combination of values of the 2 rasters. ArcGIS internally controls the resampling, and the resulting raster is 56 m in resolution.

Alternatively, one can split the 56m raster into two (getting two 28m rasters as result), resample the 28m rasters into 30 m rasters and then combine with the other 30 m raster.

Which way is better in terms of being more accurate? Any citations to this extent would be super helpful.


texts: How to save colour (color) ramps in QGIS

texts:
This strikes me that it should be simple, but efforts in QGIS and research here and elsewhere are not providing a solution.

I have a (vector) map which I have to mark out in colours according to the values of the individual sectors.  To aid the maps use for those with colour blindness issues I am only using one colour - blue.  The colours in the standard colour ramp are too close so I have made my own thus: -

New Colour Ramp
Gradient - OK
Colour 1 - white
Colour 2 - blue
Name the Colour Ramp - Test - OK

I have now got 5 colours in front of me.  Only the White is acceptable, so I double click on each colour in turn and change their values to make them more widely spaced - OK - OK.

When I print out the map all is well, but when I save that project, then open another and try to apply the colour ramp the colours have gone back to the default values in that colour ramp.  I have messed about in the style manager with no success.  Will I have to manually set all of the values for this colour ramp in each project that I wish to apply this colour ramp to?  Surely not!  There must be a simpler way.


texts: Publishing EPSG:3857 raster in GeoServer

texts:
I have two sets of services, the first one is taking a satellite product and generating a GeoTIFF to be published in GeoServer via the second service.

The code for the first service is written in python and uses gdal bindings. The GDAL version is 3.1.0 and the code runs in a docker container which image is largely based on osgeo/gdal:ubuntu-small-latest.

When I try to publish the GeoTIFF in GeoServer I get an exception.

GeoTiffMetadata2CRSAdapter::createCoordinateSystem:Only Geographic & Projected Systems are supported.

After some investigation it seems to be down to the CRS it doesn't recognize. If I use 4326, then I have no problem publishing, if I use 3857 it fails.
If I gdalinfo the faulty GeoTIFF file this is what I get:

PROJCRS["WGS 84 / Pseudo-Mercator",
BASEGEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]],
CONVERSION["Popular Visualisation Pseudo-Mercator",
    METHOD["Popular Visualisation Pseudo Mercator",
        ID["EPSG",1024]],
    PARAMETER["Latitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8801]],
    PARAMETER["Longitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["False easting",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]]],
CS[Cartesian,2],
    AXIS["easting (X)",east,
        ORDER[1],
        LENGTHUNIT["metre",1]],
    AXIS["northing (Y)",north,
        ORDER[2],
        LENGTHUNIT["metre",1]],
USAGE[
    SCOPE["unknown"],
    AREA["World - 85°S to 85°N"],
    BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]

I am using a GeoServer docker container with GeoServer 2.17, image is geosolutionsit/geoserver:cloudsi-2.17.x

How can I make the GeoServer publish my 3857 projected GeoTIFF? Do I need to add this "new" projection?


texts: ogr2ogr slow to load ArcGIS file geodatabase data into PostgreSQL

texts:
Problem: Using ogr2ogr to transfer feature sets and attribute tables from ArcGIS file geodatabases to a PostgreSQL database is taking unnecessarily long. For instance, regardless of machine I am working on, it takes roughly two days to load a table with 9 million rows; each row has five attributes with formats of double precision, integer, integer, varchar(55), and double precision. I am looking for help because this is going to be an on-going task and I am trying to develop a repeatable work process.

What I am hoping to find out: I have done some searches on speeding up ogr2ogr and PostgreSQL, but the most common suggestions have not sped things up. My guess is there is either a flag I am missing in ogr2ogr or a setting I am overlooking in postgresql.conf that I need.

What I have tried:
The ogr2ogr command format I am using is:

ogr2ogr -progress -overwrite -skipfailures --config PG_USE_COPY YES -f "PostgreSQL" PG:"host=*hostip* user=*username* dbname=*databasename* password=*userpassword*" "C:\pathTo\GDBName.gdb" "*tablename*"

Based on PostgreSQL tuning recommendations I found, I made the following changes to postgresql.conf:

shared_buffers - increased to 500MB
temp_buffers - increased to 32MB
work_mem - increased to 32MB
maintenance_work_mem - enabled but left default value (64MB)
wal_buffers - increased to 8MB
checkpoint_timeout - increased to 15min
max_wal_size - increased to 15GB
checkpoint_completion_target - increased to 0.85

Environment:

PostgreSQL computer - Windows 10 Enterprise (64-bit), Intel Xeon CPU (2.6GHZ), 16GB RAM, PostgreSQL 9.5 with PostGIS. This is not an ESRI Enterprise Database (SDE).
Data source computers (various) - So far I have used Windows 7, Windows 10, and Ubuntu 17.04 computers to run ogr2ogr. All have 8GB RAM and CPUs that are 2-3 years old. All are using GDAL 2.2.1.

**Note: I know the PostgreSQL database can load feature sets and attribute tables from a file geodatabase quickly. I first loaded my data using ArcCatalog 10.5 and everything loaded quickly. The large tables took 2-4 hours on average as opposed to the 2 days to 1.5 weeks ogr2ogr is taking. The problem with ArcCatalog, and the reason I am using ogr2ogr, is that it converted all the double precision ID fields to integer, which made the ID fields useless.

While corrupted, the fact that the database did load large (to me) tables rapidly using ArcCatalog is why I am assuming I have overlooked a simple configuration tweak.


texts: Migrating from FME Server 2008 (with Spatial Direct) to FME Server 2011?

texts:
I am about to migrate a web app from it's old architecture to new.  Part of this is migrating the old FME Server 2008 (using a Spatial Direct service) instance to FME Server 2011.

I understand it is quite different; has anyone got any experience of moving to 2011 and can offer some advice?


texts: Securing AGS output, cache, and jobs folders

texts: ArcGIS Server (for .NET) generates images, etc. in arcgisoutput, arcgiscache, and arcgisjobs folders. Can the AGS security system (specifically, token-based) be used to secure these folders in the same manner as the REST API? Or do they require a separate strategy? I've pored over the documentation but have found nothing specific to this aspect of AGS security.


texts: Updating Fields based on Selected Record/Spatial Location

texts: I have two SDE Feature classes, a point layer for Buildings and a Polygon layer for Building Footprints.  I am creating the Polygons on a job by job basis and would like to have the attributes for the point, inside said polygon, to carry over to the attribute fields of the polygon.  I am using ArcGis 10.2.  I have looked at Spatial Joins and tried field calculations, but I can't figure out how to update a record for the existing SDE layer.  I would like model to share with my team.


texts:
Is there a way to create color separated PDFs through QGIS's print composer?

I've also tried different PDF printers, but I haven't find any free option for it.


texts:
I have install QGIS 2.99 version but I have dll error missed every time to try start it.
I try different ways to do that install standalone or via OSGEO but I have some error.

I use windows 10 64bit.

error message dll miss:

qgis_app.dll
qgis_analysis.dll
qgis_qui.dll

how to fix that ?is the bug of QGIS ?


texts:
I have a scanned map of the US in pdf that, in addition to state boundaries, has "delivery territories." These territories are exclusive areas within the US where different companies deliver their products. The boundaries of these territories may cross state lines (covering at least part of two or more states), or they may be contained within a state. These boundaries are defined at the zip code level.

I have already georeferenced the pdf map and if I overlap a shapefile of zip codes, the boundaries of the territories match the boundaries of the zip codes (though territories cover many zip codes, and often a complete state or more). What I want to do is to recover the information about delivery territories from this map. Is there a way to do this? In the map, state boundaries are represented using red lines, while delivery territories are represented using black lines. I can't post  the map online because of copyright issues, but I got the following map online that represents a similar setting. In this setting, I would like to recover the territories, for example, at the zip code level, that correspond to each color (in my case, a different color would represent a different company)

Edit: My map is a white map with only the boundaries in color (red and black). The map below would be very similar to mine if instead of colors representing different areas, it only had the boundaries.


texts: Looking for an Android Object-relational mapping programming layer, PostGIS oriented

texts:
I'm working on a Android app which should manage data on a remote PostGIS database. Before diving into actual programming, I'm searching for ORMs for Android. I've found a couple, namely androrm and ormlite.

However, these don't seem to give me an interface to query objects by the geometry column. For example, find all objects within a radius from the current location... I'd need to write SQL for this...

What are my best options here?

Is there a way to completely avoid PostGIS SQL?

Or not really?

Thanks in advance.

(EDIT, 19/08/11: changed the title)


texts: Is there any way to load pre-rendered svg from the server into the Polymaps?

texts: I need to render 30000 grid cells (polygon) into a Polymap map. Is there any way to load pre-rendered svg from the server into the Polymap?


texts:
i am trying to set code for ODBC connection in Mapbasic. Instead of using Username and password i wanted to set code for Window authentication.my sample code is

Register table "tablename"
Type ODBC
Table "Select * From Database.tablename"
CONNECTION
"DSN=Database; UID=Username; PWD=Password"
into "C:\temp.table.tab"

so in this case instead of UID and PWD i wanted to set as window authentication to connect to database...anybody knows this.


texts: How to change color of a feature in openlayers?

texts:
I'm loading a geojson file into openlayers to show some polygons. From this file I also create a list of all names of these polygons.

Now I want to change the color of a polygon when selected from the list (clicked on the name).

What I've tried is to create a style (http://docs.openlayers.org/library/feature_styling.html) but I couldn't find out how to add this style to the polygon. How can I do that?


texts: Losing resolution of my DEM/TIFF when converting from HGT

texts:
I need to get a nice hillshade from this area. My steps:

1) Get HGT files from the area:
 -- phyghtmap --download-only --a -43.7544:-23.2363:-42.0378:-22.3183 --earthdata-user=myuser --earthdata-password=mypass --viewfinder-mask=1

2) Create a VRT file to join all downloaded HGT:
  -- gdalbuildvrt ./test.vrt hgt/SRTM1v3.0/S23W043.hgt  hgt/SRTM1v3.0/S23W044.hgt  hgt/SRTM1v3.0/S24W043.hgt  hgt/SRTM1v3.0/S24W044.hgt

3) Create the TIF:
  -- gdaldem hillshade test.vrt test.tif -z 10 -s 90000

The result is very ugly:

What I'm doing wrong?

EDIT

This is my GDALINFO result:

Driver: GTiff/GeoTIFF
Files: test.tif
Size is 7201, 7201
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 = (-44.000138888888891,-21.999861111111109)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -44.0001389, -21.9998611) ( 44d 0' 0.50"W, 21d59'59.50"S)
Lower Left  ( -44.0001389, -24.0001389) ( 44d 0' 0.50"W, 24d 0' 0.50"S)
Upper Right ( -41.9998611, -21.9998611) ( 41d59'59.50"W, 21d59'59.50"S)
Lower Right ( -41.9998611, -24.0001389) ( 41d59'59.50"W, 24d 0' 0.50"S)
Center      ( -43.0000000, -23.0000000) ( 43d 0' 0.00"W, 23d 0' 0.00"S)
Band 1 Block=7201x1 Type=Int16, ColorInterp=Gray
    Computed Min/Max=-38.000,2277.000
  NoData Value=-32768

EDIT 2: Using gdaldem hillshade test.vrt test.tif -z 1  -s 80000

... and this is my best result so far. After load the tif as a Geoserver layer. You can see the very poor "pixel as big sqares" result:

...same image from https://www.opencyclemap.org/ ( Pão de Açúcar, Rio de Janeiro, RJ, Brasil):

EDIT 3: Making progress:

gdal_translate -tr 0.00009 -0.00009 -r cubicspline -of GTiff test.vrt test.tif

gdaldem hillshade -co TILED=YES -co compress=lzw -s 111120 -z 1 -az 315 -alt 60 -compute_edges test.tif final.tif


texts:
using Web Feature service and I'm wondering is there a way using get requests to pull multiple features. For example, this WFS http://demo.mapserver.org/cgi-bin/wfs?SERVICE=WFS&VERSION=1.0.0 contains 2 layers/features:

Cities

Continents

Could anyone point me in the right direction for this?


texts:
I'm having problems with visualizing a layer which is in EPSG:4326

Apparently the geometry looks perfect, but I don't know why, it's not in real coordinates.

I have projected it to 4326 with QGIS, re-projected with ogr2ogr, and the result is always the same, the coordinates shown in the box while I pass the mouse are not real.

SRC of the project is in 4326, and 'on the fly' is set to on.

Any idea what could be happening?
I attatch the europe roads layer, if anybody wants to open it and check for some issues.

Thanks a lot!

Shapefile: https://dl.dropboxusercontent.com/u/64621273/roads.rar


texts: How do you get a color topographic map in Google Earth using the altitude provided by the 3D models

texts:
Is there any way to get an accurate, high definition, color topographical map as an overlay in Google Earth?  It would be somewhat similar to this image captured by NASA of the surface of the moon

Except as an overlay and of the Earth of course.  I am looking for a way to possibly generate that map from the 3d model of the earth that Google provides as I only need it for a small area and it needs to be high accuracy.


texts: Running Python produces 'No module named 'PyQt4' in Qgis on Mac High Sierra

texts:
I was following the instructions in this guide: Create a Polygon Shapefile from Point data and fairly sure I substituted field variables to match my data. However, when I run the script on my Mac (High Sierra) I get the following messages:

exec(open('/var/folders/yw/x44kj4x96mx70gych38mvzcm0000gn/T/tmpr48c2xuz.py'.encode('utf-8')).read())
Traceback (most recent call last):
  File "/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 2, in <module>
  File "/Applications/QGIS 3.app/Contents/MacOS/../Resources/python/qgis/utils.py", line 664, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)
ModuleNotFoundError: No module named 'PyQt4'

I even tried Installing PyQt4 On Mac OSX and Add to python path mac os x
but to no avail, still get the same error message.

Any thoughts here?

Here is the code I used:

import math
from PyQt4.QtCore import QVariant

def rotate(point, angle):
    angle *= math.pi / 180.
    px, py = point
    qx = math.cos(angle) * px - math.sin(angle) * py
    qy = math.sin(angle) * px + math.cos(angle) * py
    return qx, qy
def translate(point, translation):
    return point[0] + translation[0], point[1] + translation[1]

def rectangle(width, height, angle, center):
    p1 = (-width/2, -height/2)
    p2 = (-width/2, +height/2)
    p3 = (+width/2, +height/2)
    p4 = (+width/2, -height/2)
    rotated = [rotate(p, angle) for p in [p1, p2, p3, p4]]
    translated = [translate(p, center) for p in rotated]
    points = [QgsPoint(x, y) for x, y in translated]
    return points

in_layer = QgsMapLayerRegistry.instance().mapLayers().values()[0]

layer = QgsVectorLayer('Polygon?crs=epsg:2805', 'poly' , "memory")
pr = layer.dataProvider()
pr.addAttributes([QgsField("rate", QVariant.Double)])
layer.updateFields()
for feature in in_layer.getFeatures():
    geom = feature.geometry()
    p = geom.asPoint()

    width = 0.8*0.3048 * feature.attribute('width') / 100000 # MAGIC
    length = 1.2*0.3048 * feature.attribute('length') / 100000 # MAGIC
    direction = feature.attribute('angle')
    rate = feature.attribute('id_spaces')
    center = (p.x(), p.y())

    poly = QgsFeature()
    points = rectangle(width, length, -direction, center)
    poly.setGeometry(QgsGeometry.fromPolygon([points]))
    poly.setAttributes([rate])
    pr.addFeatures([poly])
layer.updateExtents()
QgsMapLayerRegistry.instance().addMapLayers([layer])

where I substituted 'width','length','angle','id_spaces' from my own point shapefile.


texts: Add additional markers to existing featuregroup in Leaflet

texts: It is possible to add additional markers to an existing featuregroup after it has already been added to the map?


texts:
I'm using GDAL v2.14 on my Windows 10 machine. In that environment, I have two gdal python bindings. One is import gdal, the other is import osgeo. I installed Python's gdal using http://www.lfd.uci.edu/~gohlke/pythonlibs/#gdal.

On my Ubuntu server, there are two ways to install gdal. One is to do sudo apt-get install python-gdal and the other is pip install gdal (but you need the C headers for this install).

On my Windows install, I can use gdal.Warp() to warp a projection, or crop it. On Ubuntu, this fails with the error AttributeError: 'module' object has no attribute 'Warp' on both "installations" of gdal. How can I make it cross-platform while also maintaining the same functionality?


texts: Error in which.max() : not yet implemented for large objects in R

texts:
I am going to create raster by Julian day that observed maximum value of NDVI within a single year. Actually, I want to know which julian day (or raster) has max value in each cells and get result by a single layer of raster.
I am working with huge raster stack (3153, 8597, 27106341, 11  (nrow, ncol, ncell, nlayers)).
I have found following script to do that, but it reported error message as like "Error in which.max() : not yet implemented for large objects".

Is there any solution for this kind of error?

# Required packages
library(raster)
library(rgdal)

# List files of interest, e.g. MODIS Terra NDVI files for 2000
fls.ndvi <- list.files("D:\\NDVI_MODIS_500\\Mod_ndvi_500_real_16",
                      pattern = "mod13a1.a2000.*.*.img$",
                      full.names = TRUE, recursive = TRUE)
fls.ndvi

# Extract julian day from available filenames
jdn.ndvi <- as.numeric(substr(basename(fls.ndvi), 14, 16))
jdn.ndvi

# Import files into R
rst.ndvi <- stack(fls.ndvi)
rst.ndvi

# Per pixel, identify layer holding the maximum EVI value
rst.ndvi.max <- which.max(rst.ndvi)
### Error in which.max(rst.ndvi) : not yet implemented for large objects

# Replace no. of layer with corresponding Julian day
rst.ndvi.max[] <- jdn.ndvi[getValues(rst.ndvi.max)]

What should I do to avoid above errors?


texts:
I'm using MakeQueryLayer to create a view from a database table in the ArcMap python window. The command completes successfully but I can't see the result. Where does it go?

arcpy.MakeQueryLayer_management(r'Database Connections\Connection to database.sde',
'csw_test','select * from REPORT_TBL',"OBJECTID (Long Integer, Not nullable)")
<Result 'csw_test'>

Nothing shows up in table of contents, List by Source is active, changing panes and forcing a refresh has no effect. Using AddTableView doesn't work:

arcpy.mapping.AddTableView('Layers','csw_test')
Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\arcgis\desktop10.6\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\arcgis\desktop10.6\arcpy\arcpy\mapping.py", line 113, in AddTableView
    assert isinstance(data_frame, DataFrame)
AssertionError


texts: Summing data from column in the Attribute table

texts:
I'm working with population data, and I have an issue when it comes to knowing the sum of a population inside the attribute table.

In my attribute table, I'd like to add all the numbers in the red box, i.e. all the numbers together from the column "UN_2015_E".
Is there a way to do so?


texts:
I have a points layer and a polyline layer.

The points layer is basically bus stops and the line layer is the road network.
The bus stops may or may not fall exactly on the road network layer.
I want to create a column in the point layer which contains the ID of the nearest link and also the distance from the link's head node or tail node.

That means I need to find the nearest link wrt each node and also measure the distance of the point on the line from the line's head or tail node.

What tool in QGIS will help me get the above things done?


texts: GDAL GetGeoTransform returning lat/lng instead of affine co-ordinates

texts:
I'm using the GDAL libraries via the c# wrapper and have hit a problem when using geotiff's where the co-ordinates are from a state plane zone. Previously I was having trouble getting the correct affine co-ordinates converted to lat/lng.

Having figured out the problem there, I now find that when I get the GeoTransform values from the dataset, it contains lat/lng rather than the affine co-ordinates. This makes it harder to calculate he co-ordinates of each corner of a raster. (In fact, I'm not sure how to do that).

On top of that, I have to somehow detect that this raster is going to give me lat/lng rather than affine co-ordinates. The API doc doesn't mention that lat/lng may turn up in the GetGeoTransform()...should it?

My code:

double[] geot = new double[6];
ds.GetGeoTransform(geot);

returns:
-77.5286920001666,
1,
0,
38.728995261289775,
0,
-1

Is there some way I can either find the affine co-ordinates using the lat/lng, or even better, get the dataset object to give me the affine co-ordinates instead of the lat/lng?

EDIT: Adding GetProjectionRef() here:

PROJCS["NAD_1983_HARN_StatePlane_Virginia_North_FIPS_4501_Feet",
GEOGCS["NAD83(HARN)",
DATUM["NAD83_High_Accuracy_Reference_Network",
SPHEROID["GRS 1980",6378137,298.2572221010002,
AUTHORITY["EPSG","7019"]],
AUTHORITY["EPSG","6152"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4152"]],
PROJECTION["Lambert_Conformal_Conic_2SP"],
PARAMETER["standard_parallel_1",39.2],
PARAMETER["standard_parallel_2",38.03333333333333],
PARAMETER["latitude_of_origin",37.66666666666666],
PARAMETER["central_meridian",-78.5],
PARAMETER["false_easting",37673535.76388889],
PARAMETER["false_northing",21527734.72222222],
UNIT["us_survey_feet",0.3048006096012192],
AUTHORITY["EPSG","2853"]]


texts: We are trying to build out stack on top of geo-spatial data, this means storing more than 10000records per a second. we are using Storm to accelerate receiving the data, but we would like to use some powerful tools like Cassandra to store data. The problem with Cassandra is that it lucks some GIS features. Do you recommend any other Data storage platform to easily query this massive data later ?


texts:
The LA Times has a leaflet map viewer app showing all 30,699 shots in Kobe Bryant's career.  It seems like spatial analysis could be useful for certain coaching decisions.  Are any coaches actually using GIS tools?

Update2

LA Times wrote an article on how they did it, with more details on GitHub.

Update
I see the NY Times has published some even more impressive analytical tools, based on work by Prof Kirk Goldsberry.

Also, BasketBallreference.com has heatmaps, like this one for Kobe.

I wonder if coaches try to get in the zone the way players do.  If so, would tools like these be a distraction?

Or does the color phi phenomenon prove that flow is an illusion?

SportVu collects the data. I don't think they were around when Kobe began his career.


texts: Mapserver GetLegendGraphic via expression

texts:
I am trying around this for a while now. But can't find any answer. So I try it here.

What do I have:

Mapserver 7.0.6

What do I want:

Obtaining a single legend graphic for a class without text and in custom width/heigth by passing the value evaluated in the class.

Let me describe it a little more with examples:

Lets say we have a CLASS defined in a mapfile with the NAME "Test Name" which matches the EXPRESSION "100". I can access this simply by the following URL:

test.mapserv?width=72&height=36&layer=test&version=1.1.1&service=WMS&format=image/png&request=GetLegendGraphic&rule=Test Name

I will get the resulting image as I expect. But for the reason of encoding issues and because I trust the numbers more than prosa text I want to use the expression instead. My dream would be something like this:

test.mapserv?width=72&height=36&layer=test&version=1.1.1&service=WMS&format=image/png&request=GetLegendGraphic&rule=100

For sure: This won't work. I tried several other stuff. I was passing a SLD to filter the legend by the EXPRESSION but this answers with the legendimage + the legend text like the following:

This is explicitly not what I want. I only want the simple icon. And I really want to do it via URL of OGC service to stay as close as possible to well known standards.

So now my hope is to find somebody out there who is able to show me a way to achieve this.


texts: ArcGIS network analyst load locations is very slow

texts:
I have been using a location allocation model for several months now, with no issue. The load location feature is extremely slow, minutes for 1 point.

There is a spatial index and i have deleted and updated them several times. I have tried older versions and used older (faster) networks. All to no success.

I have seen this question and answer:

Load locations faster in ArcGis network analyst

Is there anything else to be done?


texts: ArcMap 10.2.1 issues diffferent SQL than 10.0 for "Select by location"?

texts:
I have a large point layer and I'm doing a "Select by Location" to select all points that intersect a specific selected polygon from a polygon table.
In ArcMap 10.0 the SQL that was executed was basically the following:

Stick the se_rowid of the polygon into the sde_logfile_data table

Select all points that intersect the polygon

In ArcMap 10.2 the SQL is an alternating sequence of:

Insert the rowids of 10000 points into the sde_logfile_data table

Select count(distinct sde_row_id) from sde_log_table where logfile_data_id = ...
i.e., sde_logfile_data keeps getting inserted into and the subsequent selects take longer and longer.  ArcMap ends up stalling and I have to kill it.

Does this sound familiar to anyone?


texts: How can I create multipolygon using shapely?

texts:
I am going to create polygon, but it produce invalid polygon form this data:
[[  0 438]
 [  0 428]
 [  0 427]
 [  6 428]
 [  9 434]
 [ 11 439]
 [  7 443]
 [  3 446]
 [  0 448]
 [  0 452]
 [  0 452]
 [  0 449]]

Now I want to create multipolygon from the above data. How can I create multipolygon?


texts:
I want to associate a point layer with a polygon layer so that I can do zonal statistics. In the image below, the polygon layer represents roughly the metro area of a certain city. The point represents the city. The polygon is a result of a polygonized process from a raster image. I want to use these new polygons as buffers to analyze the original raster. However, I need a way to associate the city points with the surrounding polygons.

In other words, I want to find zonal statistics for each city ("oecdminusa") using the corresponding metro polyon ("polygontest"). When I ran zonal stats the first time, it simply named each of the polygons in "polygontest" as 1,2,3,4 etc. How do I make it correspond with the appropriate city?


texts: Writing domains with coded values in ArcSDE Enterprise Geodatabase

texts:
I'm trying to write a new domain and values in our ArcSDE database but I receive the following error:

Runtime error  Traceback (most recent call last):   File "", line 34, in    File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\management.py", line 1317, in AddCodedValueToDomain     raise e ExecuteError: ERROR 000800: The value is not a member of Gas Service Diameter Plastic PVC | Rectifier Manufacturer | Work Order Status | GasASTMF2897CSteelPipeType | Pipe Diameter Plastic PVC | Leak Survey Frequency..

Line 34) arcpy.AddCodedValueToDomain_management(gdb, domName, code, domDict[code])

Notes: The domain name is successfully created in the enterprise geodatabase but the coded values return an error.

Current ArcSDE Enterprise Geodatabase Version: 10.3.1.

Code:

#Import system modules
import arcpy

arcpy.env.workspace = "C:/Users/jnmiller/AppData/Roaming/ESRI/Desktop10.3/ArcCatalog"

# Set local parameters
domName = "test123"
gdb = "GIS.sde"

# Process: Create the coded value domain
arcpy.CreateDomain_management(gdb, domName, "Description", "TEXT", "CODED")

# List of domain values in load order
domainList = ["1"]

# Store all the domain values in a dictionary with the domain code as the "key" and the
# domain description as the "value" (domDict[code])

domDict = {"1":"test 1"}

# Process: Add valid material types to the domain
# use a for loop to cycle through all the domain codes in the dictionary
for code in domainList:
    arcpy.AddCodedValueToDomain_management(gdb, domName, code, domDict[code])

print "Finished!"


texts:
I'm new to PyQGIS and I'm having some dificulties. I'm trying to implement / utilize "Add Feature" (polyline) functionality in my python plugin.
I know I can do this by defining starting and ending points.

PyQGIS, construct a polyline:

line_start = QgsPoint(50,50)
line_end = QgsPoint(100,150)
line = QgsGeometry.fromPolyline([line_start,line_end])

This is what I'm trying to avoid. I'd like user to digitize / create polyline by clicking on map canvas to preview geometry. How can I do that?


texts: QGIS and MapInfo Projections - EPSG 27700 (British National Grid)

texts:
I have some data in a PostGIS table with the projection EPSG 27700. If I export this data from QGIS (v1.8) to MapInfo TAB / MIF when I open the data in MapInfo it has a shift on it.

Does anyone know why this happens? I don't get any problems if I open MapInfo data directly into QGIS, but I did notice that QGIS opens the data with a custom projection.


texts:
I need to find a way to merge ways into one long lines with the same attributes. I heard I can use ST_Union but I don't know how. I have a postgresql database connect with tilemill. It works fine but I'm not familiar with postgis so I don't know how to write the subquery. What I want to achieve is to merge segments of ways into one.

Here is my subquery in tilemill:

( SELECT *
    FROM planet_osm_line
    WHERE route IN ('bicycle')

) AS data

And here is my table in postgresql:

These are ways converted from bicycle routes. I want to combine them into one. The reason is that when rendering every single segment gets split and on angles it doesn't look good.

http://pomorskieszlaki.pl/#16/53.7062/16.6908/osm-rowerowe


texts: Combine / Group - ArcGIS attributes table

texts:
I need to group several rows in a single one, combining data cell (mainly “text”) in Arcgis desktop 10.2.2. One cell could have several records separated by comma. I have attached a screenshot in order to explain what I want to do. I have used geoprocessing tools, but I can’t figure out this problem.

I thought that one solution could work with table in excel, but I didn’t know how to make it yet (not manually).


texts: Save millions of points in the server side

texts:
I have asked a similar question about millions of points a few days ago:Display millions of points in the browser.

It focused on how to display these data in the client(browser), however I think it is still a big challenge to save the in the server side.

For example for a device like a truck, which will send its location to the server every 60 seconds(not short for trucks), and suppose the truck just work 15 hours one day, so there will be 15 * 60 = 900 records one day.

And if there are 100 trucks in the system, there will be 900 * 100 * 30 = 2700000 = 2.7 million points one month.  Though we can reduce the points by removing the repeat points, but this can not take a explicit effect.

Now we use mongodb to save these points, however the storage size increase rapidly.

I wonder how do you handle this kind of use case?


texts: PostGIS - ST_Within or ST_Disjoint performance issues

texts:
I have a table of about 150,000 points in PostGIS with a spatial index and a SRID of 27700 (OS).  I want to select the points that fall outside England and Wales.  I have a multipolygon table with just one record in it (a dissolved polygon of England and Wales).

I would have thought that the following query would select this quite quickly, but it was still running after about 5 hours (I actually forgot I'd started it and was off doing something else).  I've tried a few variations using ST_Within to do it the other way around, but I can't get a result in a reasonable time (I have yet to wait for the query to finish running).  Considering the points were generating in pgRouting in about 5 mins, I'm surprised it takes this long to do any analysis of them.

Admittedly I'm using an old PC but still with 2GB of ram I would have thought this was a relatively simple query.  Is it likely that it should take this long or is something likely to be wrong. The query I've used is below:

SELECT pk
FROM catchment_distance_output, england_wales_os
WHERE ST_Disjoint (catchment_distance_output.geometry, england_wales_os.geometry);

I imagine I'm doing something obvious, but I'm very new to this and learning as I go.


texts: Merging layers in cost distance analysis to account for other variables?

texts:
For example, say I was rafting down a river but there was a small bridge in the way. That would significantly change the difficulty but not necessarily stop the progress.

Is there a way for me to include polygon layers (from a buffer) or lines or points into the cost path to represent a change in terrain?

I'm using ArcMap 10.6.1.


texts: Qgis not computing raster stats in layer properties - metadata?

texts:
I used to check the sum of all values of a raster (to get some volume information) from the layer properties - Metadata - properties tab, where I had directly the figure I was looking for.

However, since a few weeks, I realised that the layer statistics are not computed any more.

Is there a way to tell Qgis to compute these data?

I have the feeling that this is a bug linked with the 2.18 upgrade, as I realised the misbehaviour a short time after it. Is this general or just happening on my side?

It might also be worth mentioning that I removed a few plug-ins in the upgrade process. I've not seen any plug-in computing these values in the catalogue and would assume it is a core function, but we never know.


texts: Error 000732 - input raster does not exist when implementing script tool arcpy

texts:
I am trying to implement a script tool in ArcGIS 10 that completes a few equations with raster datasets.

I seem to be getting the same error no matter what I try:

<type 'exceptions.RuntimeError'>: ERROR 000732: Input Raster: Dataset 5 does not exist or is not supported
Failed to execute (TWI).

TWI being the name of the script tool.

My code is below:

import arcpy
from arcpy import env
import os
from arcpy.sa import *
arcpy.CheckOutExtension("spatial")

DEM = arcpy.GetParameterAsText(0)
cellsize = arcpy.GetParameterAsText(1)
FlowAcc = arcpy.GetParameterAsText(2)
K_Sat = arcpy.GetParameterAsText(3)
OutRaster = arcpy.GetParameterAsText(4)
OutSpace = arcpy.GetParameterAsText(5)

slope = (Slope(Raster(DEM)) * 1.570796 ) / 90
tan_slp = Con(slope > 0, Tan(slope), 0.001)
FlowAcc_scaled = (Raster(FlowAcc) + 1 ) * cellsize
TWI = Ln(FlowAcc_scaled / (Raster(K_Sat) * tan_slp))

arcpy.env.workspace = OutSpace

TWI.save(OutRaster)

I've got a feeling that the error is something to do with setting up the parameters correctly. It's something that I seem to be struggling to get my head around but as I understand it "DEM", "FlowAcc" and "K_Sat" are raster layers, "cellsize" is a long integer (also tried double), "OutSpace" is a workspace and "OutRaster" is a string.

The error seems to be with the way is is saving/naming the raster.

EDIT:
Please see a screenshot of the parameter setup below (note I have tried the output workspace and output raster as both inputs and outputs and I get the same error)

EDIT 2 - Answer

Ok the issue was that the second parameter 'cellsize' needs to be called with the 'arcpy.GetParameter' function rather than the 'arcpy.GetParameterAsText' function as it is an input value rather than a string.


texts: After converting .shp to .gpx, there seems to be no content in the .gpx

texts:
I tried to convert .shp files into .gpx using QGIS. Originally, the .shp were in the Swiss Coordinates (CH190) but I've changed them into WGS.

I can load the .gpx into QGIS which works really fine. Trying to load it into Movescount / Google Earth / GPStracks.com, it fails.

I tried to with

GPS_use_extensions YES and NO.

I checked

skip attribute creation

and I tried

FORCE_GPX_TRACK YES and NO.

Any ideas what I do wrong?


texts: Add multispectral raster data to a mosaic dataset

texts:
I have MS images (G, R, RedEdge, NIR) and want to add them to a mosaic dataset. When creating the mosaic dataset, I cannot find an appropriate 'product definition' for it. I tried manually adding bands and wavelengths but end up with a reflectance map with two bands.

I also tried selecting 'Landsat MSS' as product definition and followed and tried to find my way by using a tutorial with the example for Landsat 7 ETM+ but it is clearly not suitable for my case.

I hope my mission is defined clearly :)


texts:
I have four markets 1000-4000, but I only want to do markets 2000 and 3000. I created a markets list only containing those two markets so that I could use the list in a "for" loop. However, when I try and run it, it still runs through all the markets. I think this has to do with my "whereclause" or "sqlclause" or both, but I can't figure out how to fix it. I've tried using a wildcard in the "whereclause" but that doesn't work. I've also tried getting rid of the "whereclause", as well as an "AND" but that also doesn't work. What is wrong?

Code:

# Import system modules
import arcpy
import datetime
from datetime import timedelta
import time
import calendar

# Set environment settings
arcpy.env.workspace = r"C:\arcGIS_Shared\Python\CenterHeatMaps2.gdb"

#Declare variables
fc = 'Open_GoHealth_Centers'
fields = ['USER_market_id','USER_GoHealth_ID','USER_GoHealth_Center_Name','USER_Opening_Date']
fieldname = 'USER_market_id'

#Define WHERE clause statement
whereclause = """{} = 2000 AND 3000 """.format(arcpy.AddFieldDelimiters(fc, fieldname))

sqlclause = (None, 'Order By USER_market_id, USER_GoHealth_ID')

# Set Markets to loop through
markets = [2000,3000]
# Set Years to loop through
years = [2016,2017]
# Set Months to loop through
months = [1, 2]

#returns last day of each month
def last_day_of_month(any_day):
    next_month = any_day.replace(day=28) + datetime.timedelta(days=4)
    return next_month - datetime.timedelta(days=next_month.day)

for market in markets:
    print (market)
    for year in years:
        print (year)
        for month in months:
            print (month)
            with arcpy.da.SearchCursor(in_table = fc, field_names = fields, where_clause=whereclause, sql_clause=(None, 'ORDER BY USER_market_id, USER_GoHealth_ID')) as cursor:
                #Loop through each row established in cursor
                for row in (cursor):
                    # Set local variables for the FeatureClasstoFeatureClass tool
                    inFeatures = "PatientVisitsGeocoded"
                    outLocation = r"C:\arcGIS_Shared\Python\CenterHeatMaps2.gdb"
                    outFeatureClass = "PatientVisits{0}_{1}_{3}_{2}".format(row[0], row[2], month, year)
                    delimitedfield = arcpy.AddFieldDelimiters(arcpy.env.workspace,"USER_CenterID")
                    expression = """{0} = {1} AND USER_DOSMonth = '{2}' AND USER_DOSYear = {3}""".format(delimitedfield, row[1], calendar.month_abbr[month], year)

    #               Execute FeatureClassToFeatureClass tool
    #               This If statement determines if the Opening_Date of the respective center is before the last day of the month being run.
    #               If it is, the feature class is copied, if not, the feature class is skipped.
                    if last_day_of_month(datetime.date(year, month, 1)) >= row[3].date():
                        arcpy.FeatureClassToFeatureClass_conversion(inFeatures, outLocation, outFeatureClass, expression)


texts: Dissolving two different buffers together so not overlapping using ArcGIS Desktop?

texts:
I am trying to find out how many registered sex offenders live within 2,000 feet of a school/park in a specific county. The park buffer and school buffer are dissolved but not together. I am also trying to find the area in acres of the 2 combined buffer zones.

I am using ArcGIS Desktop 10.


texts: Getting geometry column integer SRID from a GeoPandas geodataframe?

texts:
The primary answer to Adding GeoPandas Dataframe to PostGIS table? requires entering the geodataframe geometry column's integer SRID.

Getting the CRS attribute using the crs attribute like so

import geopandas as gpd
file_path = 'my_geodata.shp'

gdf = gpd.read_file(file_path)
geom_srid = gdf.geometry.crs['init']

gives a string like 'epsg:32616'.

The integer value can be extracted using regex:

geom_srid_num = re.search(r'\d+$', geom_srid)

But I'm not sure if the SRID is written in this format all the time, and if this method is reliable. Is there a more appropriate way to retrieve the integer SRID for a geodataframe geometry column?


texts: let user add his location on interactive map

texts: I would like to know if it is possible to create an interactive map where the user will be allowed to add their own location to the map (I would use OSM as the basemap) and then calculate the distance between that location and some other objects that the map shows? If it is possible I would appreciate if you could explain to me the basic idea how it could be done. Thanks


texts: Check fields data type before insertion to table

texts:
I am trying to insert a python list of tuples to arcgis table. there is type mismatch which is obvious. Each Value in tuples have different type and fields in table have different data types.

So, how can I loop through each field, check data type, perform type cast on each values in python tuples and insert them in table. Or is there easier way to do this ?

   pylist= [(1.0, u'10',''), (2.0, 1.0, 2.0), (u'9', '', 2.0)]
   fields =["A", "B", "C"]  #fields of arcgis table

  row = arcpy.da.InsertCursor(tablename, fields )
  for i in range(len(pylist)):
  # here i need to have a routine which checks the data type of the field and perform type cast on tuple value before insertion
        row.insertRow(pylist[i])

Any suggestion will be helpful.


texts:
Currently, I am creating GeoJson files using QGIS and render them in Mapbox GL JS. I create 10 GeoJson files for each map. I am wondering what is the best practice to render the files in Mapbox GL JS. Is loading the GeoJson files in Mapbox GL JS is good enough?
Should I use some tile server or convert the files to mvt?


texts: Unusual polygon shapes when using tree_hulls?

texts:
The following script uses lidR to segment trees and output tree hulls (i.e. polygons representing individual tree canopies). Input data are a directory of laz files and a catalog is used to manage the processing.

library(lidR)
library(rgdal)

my_tree_detection_method <- function(cluster, ws)
{
  las <- readLAS(cluster)
  if (is.empty(las)) return(NULL)

  las_n <- lasnormalize(las, tin())
  chm  <- grid_canopy(las_n, 0.25, pitfree(c(0,2,5,10,15), c(0,1), subcircle = 0.2))
  algo = watershed(chm, th = 4)
  trees <- lastrees(las_n, algo)

  # Remove the buffers
  trees_no_buffer <- lasfilter(trees, buffer == 0)
  return(trees_no_buffer)
}

ws <- "/path/to/laz/directory"
ctg <- readLAScatalog(ws)
lidR:::catalog_laxindex(ctg)

opt_chunk_buffer(ctg) <- 10
opt_chunk_size(ctg) <- 200

opt <- list(need_buffer = TRUE)
output <- catalog_apply(ctg, my_tree_detection_method, ws = 5, .options = opt)
output <- do.call(rbind, output)

hulls  = tree_hulls(output, func = .stdmetrics)
writeOGR(hulls, dsn = "/path/to/dir", layer = "hulls", driver="ESRI Shapefile")

When I write the hulls to shapefile to inspect the results, the hulls appear very distorted. I would expect thousands of small polygons representing tree canopies. What might be causing this distortion?

The following from the documentation is what I would expect to see:


texts: QGIS 3.10 (and 3.16) Mac; bad resolution and vanishing WMS/WMTS Layers

texts:
Hello dear QGIS Community,
I am setting up QGis on my MacBook M1 with Big Sur and facinging various issues with layers and external sources. After having worked for a while with QGIS on Windows 10 and having recently switched to an M1 Mac Book Pro not supporting Bootcamp, I started to use QGis directly in macOS (11.1) and via Rosetta.
Unfortunately, the setup seems to be complicated.
The resolution of connected WMS/WMTS Layers is unusable. Most of them disappear after zooming in at more than 1:12000 or so.

It concerns public sources like Openstreetmap or Openrailwaysmap as much as non-public sources.
I tried it on different display and various versions of QGIS.
It happens on every display, whether external or internal and on 3.10 as much as 3.16.

I'd be happy if someone can assist me put on me on the right direction to fix this issue - or to find out whether this is a bug at the moment due to M1 Silicon platforms.


texts: How do I export a KML for one layer in a KML file that has multiple layers?

texts:
I'm accessing the Index plots for Ordnance Survey map sheets, provided as a KML by the Charles Close Society here: https://www.charlesclosesociety.org/indexes (scroll to the bottom of the page for the KML download).

And loading the file to Google Earth, where it is possible to view all the Index plots overlayed in one big jigsaw, or view one at a time.

I would like to export one of the layers only, import to ArcGIS 10.+ and convert from KML. Although I'm viewing only one layer when I try to export, I'm still getting all the data in the exported file. Rather than the single layer I need.


texts: Hot spot analysis widget application programmed in JavaScript API is having trouble connecting to my portal

texts:
//Define the default inputs for the Hot Spot widget
        var hotSpotParams = {
            id: "analysisTool",
            analysisLayer: sex_offenders,
            portalUrl: "http://utdepps.maps.arcgis.com/home/organization.html",
            aggregationPolygonLayers: [sex_offenders],
            boundingPolygonLayers: [sex_offenders],
            showHelp: false,
            showSelectAnalysisLayer: false,
            showCredits: false,
            map: map,
            returnFeatureCollection: true
        };

        var hotSpots = new FindHotSpots(hotSpotParams, "hotSpotDiv");
        hotSpots.startup();

        //If any errors occur reset the widget (Not Working...troubleshoot)
        on(hotSpots, "job-fail", function(params){
            resetTool("error");
        });
        on(hotSpots, "job-status", function(status){
          if(status.jobStatus === 'esriJobFailed'){
          alert("Job Failed: " + status.messages[0].description);
            resetTool("error");
          }

        });
        on(hotSpots, "job-cancel", function(){
            resetTool("error");
        });
        on(hotSpots, "job-submit", function(result){
            //display the loading icon
            domUtils.show(dom.byId("loader"));
            resetTool("submit");

        });
        //The hot spots analysis has finished successfully - display the results
        on(hotSpots, "job-result", function(result){

            //hide the loading icon
            domUtils.hide(dom.byId("loader"));

            //add the results to the map and display the legend.
            if(result.value){
                var template = new InfoTemplate("Results", "${*}");
                var resultLayer = new FeatureLayer(result.value.url || result.value, {
                    infoTemplate: template,
                    outFields: ["*"],
                    opacity: 0.7, //default too transparent
                    id: "resultLayer"
                });

                map.addLayer(resultLayer);
                //refresh and display the legend
                domUtils.show(dom.byId("legendContainer"));
                legend.refresh([{layer: resultLayer}]);
            }
            if(result.analysisReport){
                //hide the hot spots panel and show the analysis info.
                domUtils.hide(dom.byId("toolPanel"));
                domUtils.show(dom.byId("infoPanel"));
                dom.byId("analysisInfo").innerHTML = result.analysisReport;
            }

            //re-enable the hot spots tool
            hotSpots.set("disableRunAnalysis", false);

        });
    }

    function resetTool(status){
        if(status === "error"){
            var tool = registry.byId("analysisTool");
            tool.set("disableRunAnalysis", false);
        }
        //Analysis started lets clean up any old results
        //Hide the legend, remove the results layer
        var layer = map.getLayer("resultLayer");
        if(layer){
            map.removeLayer(layer);
        }
        domUtils.hide(dom.byId("legendContainer"));
    }
    function showToolPanel(){
      //expand the right panel to display the content
        var cp = registry.byId("rightPanel");
        domStyle.set(cp.domNode, {width: "20%"});
        registry.byId("mainWin").resize();

        map.reposition();
        map.resize();
    }
  });
});

I'm trying to use JavaScript API to program a widget to give generate a hot spot analysis of sex offenders. My program seems to be having trouble connecting to arcgis.com (if I use https://www.arcgis.com, it prompts for a login screen. No credentials will work). I can load the feature classes, but I can't seem to get the analysis function to work.


texts: r.statistics throws invalid reclass rule error despite assigning null to unique value

texts:
Trying to run r.statistics on two raster layers, one a layer of distinct zones produced via r.clump (stands) and one a layer of all ones (ones). The areas outside of the immediate region of interest are NULL values in the original data.

Below are the lines of code:

r.in.gdal input=ones.tif output=ones -o
r.in.gdal input=stands.tif output=stands -o

r.null map=ones null=999999999
r.null map=stands null=999999999

r.statistics base=stands cover=ones method=sum output=forest

However, r.statistics keeps producing this error:

WARNING: Can't have null on the left-hand side of the rule
ERROR: 140725237898080 = 140725237898080 0 - invalid reclass rule

How do I set the nulls such that I can run r.statistics on the two layers?


texts: Hide meridians and parallel lines from plot and add TICK marks in matplotlib?

texts:
I want to hide the meridian and parallel lines from my plot and add the tick marks.
I could remove the lines by setting linewidth=0, but how to add the tick marks with each Lat-Lon label.

meridians = np.arange(70.,90.,0.5)
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=9, linewidth=0)