texts: Reclassifying a polygon using ArcGIS field calculator?

texts:
I am trying to reclassify a polygon into three classes:

1 and 2 = 5
3 and 4 = 3
5 = 1

I have been trying to do this by calculating a new field in the attribute table but I am very inexperienced with Python so have tried it VB script. Is there anyone out there that could help me with a basic IF function?

I have tried this:

if [Sal_AlkCL]=1,2 Then
 [Salscore] = "5"
if [Sal_AlkCL]=3,4 Then
 [Salscore] = "3"
if [Sal_AlkCL]=1 Then
 [Salscore] = "1"
end if


texts:
I'm trying to filter esri's geocoder for a specific country.  I've had some success, but nothing works extremely well.  I would have thought simply adding sourceCountry = "xxx" to the parameters of the arcgisGeocoder would work, but I couldn't get it to zoom correctly.   here are the three versions I've found:

attach a condition to the URL
go with no filter
use suffix = ", USA"

I'm attaching a code snippet as well below of these.

I'm wondering: Has anyone found a better way to filter the autocomplete options and have the map zoom to the correct coordinates?

//1. this one works okay, but the autofill is much slower than the others

   var myGeocoders = [{
      url: "http://geocode.arcgis.com/arcgis/rest/services/World/GeocodeServer?sourceCountry=USA",
      name: "EsriWorldLocator",
      placeholder: "Find a place",
      autoComplete: true,
     singleLineFieldName: "SingleLine"
}];

var geocoder = new Geocoder ({
    map: map,
    arcgisGeocoder: false,
    geocoders: myGeocoders
    }, dom.byId("search")
    );
geocoder.autoComplete = true;

//2. This is the original, no filter

var geocoder =  new Geocoder({
  arcgisGeocoder: {

    placeholder: "Find a place"
  },
  autoComplete: true,
  map: map
}, dom.byId("search")
);

//3.  This one has a great autocomplete, but it zooms to a seemingly random location

var geocoder =  new Geocoder({
  arcgisGeocoder: {
    //sourceCountry: "USA",
    suffix: ", USA",
    placeholder: "Find a place"
  },
  autoComplete: true,
  autoNavigate: true,
  map: map
}, dom.byId("search")
);


texts: Using Python to Union two polygons from separate gdb?

texts:
I'm in the processes of learning python and am stuck on an exercise. One part of the exercise requires me to union a polygon from one gdb with a polygon in another gdb.

Here is what I have, which isn't working:

arcpy.Union_analysis(in_features="polygonA #;polygonB #",out_feature_class="C:test.gdb/union",join_attributes="ALL",cluster_tolerance="#",gaps="GAPS")

So polygonA is from one gdb, while polygonB is from test.gdb.
I feel like I should put the pathway for polygonA in, but I'm not quite sure how to do so.

I'm using ArcMap 10.2.2.


texts: Splitting feature class by unique values for an attribute/field in ArcGIS for Desktop?

texts:
I would like to split a single shapefile into multiple shapefiles by using unique values for one attribute/field to select features for each value out into a new feature class.

I found ET Geowizard to do this but it's not free.

In arcscripts.esri.com I found a tool named "SplitLayerByAttributes" but it does not work in ArcMap . (see comments below, this was an old version which has been updated.  The new version can be found at resources.arcgis.com   )

Do you know of any tools or solution to split by attribute in the way described above?


texts: how to change encoding of a geopackage in QGIS 3

texts: For other formats you just go to Label >> Properties >> Source and there's a clickable button to select from a list of locales; but this box is grayed out for GPKG, is there a way to change the encoding of this kind of files?


texts: Remove empty spaces with FME PythonCaller

texts:
I have 4 fields in the attribute table. After I merge these attribute fields in one single attribute field using FME Desktop 2020.1, it empty records within the fields are producing empty spaces and I want to get rid of empty spaces.
what I have:

Attribute field A

1,2,3,4
1, ,3,4
 ,2, ,4
 , , ,4
 , ,3,
1,2, ,4

What I want:

Attribute field A

1,2,3,4
1,3,4
2,4
3
1,2,4


texts: Deleting polyline vertices coincident with points in ArcGIS?

texts:
I'm looking to take a shapefile of points and a geodatabase feature class of lines as input and delete the vertices from the lines coincident with the points.

I found a similar question already asked and answered, but dealing with lines and polygons instead of points. Would the method outlined here by Paul work with points as well? I have never used difference so I don't know how it would work when run on points.


texts: Is there a way to create buffer for grid coordinates in Print Composer?  I prefer to add the grid co-ords inside the map but the coords are often overwritten by map data which makes them impossible to read.  Often the coordinates are of more interest to the map reader than the data which overwrites them.  If they were buffered they would be easier to read.


texts: ArcGIS Javascript API not loading portal web map

texts:
I'm trying to set up the web map example using our portal, but nothing is actually loading

The web map is shared to public

HTML Code :

<!DOCTYPE html>

  <meta charset="utf-8">
  <meta name="viewport" content="initial-scale=1,maximum-scale=1,user-scalable=no">
  <title>Eletrosul test map</title>

  <link rel="stylesheet" href="https://js.arcgis.com/4.0/esri/css/main.css">
  <link rel="stylesheet" href="https://js.arcgis.com/4.0/esri/css/calcite/calcite.css">

Our Portal (It's in Portuguese, but the link and ID checks) :


texts: Finding optimal network with python and network analyst in arcgis 10.0

texts:
Hello GIS StackExchange,

I'd like to know if there's any way to find the shortest distance across the diagonals of a polygon. I am creating a network analyst solution to determine the most cost effective route between 2 points. The solver will travel along the roads edge for the most part to determine optimal routes but there may be situations where cutting across certain polygon features will be more efficient than continuing along the edge of the road. The question here is what is the best way to cut across the traversable polygon features?

The image is a sketch i did to illustrate how i'd like the route optimizer to function. The blue point feature is the route origin point while the green feature is the destination point. The green polygons are parcels that can be traversed, while the yellow polygons cannot be cut across. The blue line is the road network which is the main edge feature to determine the optimal route, the purple diagonal(ish) line is the result of the router determining that cutting across the polygon would indeed be more optimal than following the road.

The idea i have for cutting across the polygon features is to first determine the shortest distance within all traversable polygons and then create a line feature class containing the shortest distances for all the traversable polygons in the network dataset. This line feature class will then be connected to the same connectivity group as the roads feature class to ensure the router can connect both features and possibly determine that cutting across the polygon is indeed more optimal than following the road edge in certain situations.

One problem i have encountered in my research is that polygons may be made up of countless number of vertices and as such it would be difficult to determine the shortest path across the polygon. This question is related to my thesis project, it is to be written in python within arcgis 10.0. I have done quite a lot of research already so don't feel obliged to remind me about the need to learn this for myself because i have and i'm stuck. My biggest fear is designing the model half way only to find out that the approach i took was wrong from the start. That's why i need some help. I'd appreciate any documentation that can help or any complete resource that will prevent me re-inventing the wheel, or even ideas about how to go about certain aspects of building the router.

Thanks as always gis stackexchange. Apologies for the lengthy post.


texts: Free alternative to ArcGIS World Geocoding Service for reverse geocoding

texts: I'm seeking a free alternative to ArcGIS World Geocoding Service for running geoprocessing tool Reverse Geocode in ArcGIS Pro. All feature types provided in the ArcGIS World Geocoding Service would be ideal


texts:
I have a set of vector based shape files, representing different types of ground coverage in a city (for ex: roads, buildings, open field,...).

For the algorithm that I want to run on it, I need just one raster file with these different ground coverages represented by different bands.

The final file needs to be from IMG file type.


texts:
I am having problem when I am Requesting following URL in GeoServer

/geoserver/wfs?REQUEST=GetFeature&VERSION=1.1.0&typename=Final&propertyname=name&service=WFS&

I am getting following error:

<ows:ExceptionReport xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows" version="1.0.0" xsi:schemaLocation="http://www.opengis.net/ows http:///geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd">
<ows:Exception exceptionCode="NoApplicableCode">
<ows:ExceptionText>
java.lang.RuntimeException: Unable to obtain connection: Connection refused. Check that the hostname and port are correct and that the postmaster is accepting TCP/IP connections. Unable to obtain connection: Connection refused. Check that the hostname and port are correct and that the postmaster is accepting TCP/IP connections. Connection refused. Check that the hostname and port are correct and that the postmaster is accepting TCP/IP connections. Connection refused: connect
</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>

Can any one help in resolving this issue?


texts: Dealing with empty feature classes in a ArcGIS 9.3.1 Python script

texts:
I have a script that goes over all the feature classes within a geodatabase (in ArcGIS 9.3.1) and copies them. The problem (as far as I can understand) is that when the script goes to a  empty feature class, the script fails and gives the following error:

: Failed to execute. Parameters are not valid.
ERROR 000732: Input Features: Dataset Parcels_border does not exist or is not supported
Failed to execute (CopyFeatures).

How can I get the code to skip over empty feature classes? I tried to do a GetCount() with a if clause - if empty then skip and so on - but even the GetCount() doesn't deal well with empty feature classes either(it puts out the same error).

Does anyone know how I can overcome this error?

EDIT: Hearing to Matt's advice, here's my script:

    # Import system modules
import sys, string, os, arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(9.3)

# Load required toolboxes...
gp.AddToolbox("E:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Data Management Tools.tbx")

# Local variables...
cop_dictionary = {"xxx":"rozetta","yyy":"height_line","zzz":"walls_fences"}

part_path = "I:\\Users\\john\\"  #setting part_path

for iter in cop_dictionary:
    path = part_path + cop_dictionary[iter] + ".shp"
    if (gp.GetCount(iter)!=0):
        gp.CopyFeatures_management(iter, path, "", "0", "0", "0")
        gp .AddMessage(iter + "has been exported")
    else:
        gp .AddMessage(iter + "Is Empty")

XXX,YYY,ZZZ are feature class (FC) names in a non-english font taken from a geodatabase which is opened in the ArcMap TOC. The non-english font is not a problem, as it runs through the first FC (XXX), and exports a intact shapefile.  The error above appears in "YYY" which is a existing and empty FC.

Thanks a million.

SECOND EDIT: tried the 2 following answers below, both good (and teaching), but none of them solved the problem. It seems the problem isn't empty FC, but FC with some thing wrong in them (I can't pinpoint what). I haven't solved the problem , but left the problematic FC out until I figure out the problem.


texts:
I have exported a file from Maps.me as a KML and want to edit the style in QGIS and then export it back to Maps.me.

It is possible to edit the styleUrl in the raw KML-file, then save the changes and import it into Maps.me, but there is no option to edit the styleUrl in QGIS.

In QGIS the styleUrl field does not appear when the KML is loaded into QGIS. I was wondering if it is at all possible to edit the styleUrl in QGIS?


texts: Assigning an OSM originated tag to a Shapefile field in QGIS

texts:
I am using QGIS to do some basic data processing.

I have a shapefile (originated from OSM data) that has a number of fields. One of those fields (tags) contains a bunch of tag data.

Example of a "tags" field record below:

"name"=>"Felice Island", "place"=>"island", "osm_uid"=>"14512", "osm_user"=>"DENelson83", "way_area"=>"1.94426e-05", "osm_version"=>"1", "osm_changeset"=>"59732201", "osm_timestamp"=>"2018-06-11T08:10:52Z"

What I want to do is create a brand new field that contains just the values of the "way_area" tags.

I've been trying to set this up using the QGIS attribute table but can't find a way to do it.


texts: Hide/Show route layer in Leaflet when click on button (Easy Button plugin)

texts:
I have a Leaflet map with button (Easy Button plugin https://github.com/CliffCloud/Leaflet.EasyButton). When click on button on map added route layer (plugin https://github.com/perliedman/leaflet-routing-machine). Please, help me add event for button to show/hide layer when I click on button. Map example:

var map = L.map('map', {
    center: [50.0669, 35.1638],
    zoom: 15
});

L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 18,
    attribution: '&copy; Map Data OpenStreetMap contributors'
}).addTo(map);

L.easyButton('fa-level-up',
    function() {

        var routing = L.Routing.control({
            plan: L.Routing.plan([
                L.latLng(50.07132, 35.14103),
                L.latLng(50.05459, 35.18239)
            ], {
                waypointIcon: function(i) {
                    return new L.Icon.Label.Default({
                        labelText: String.fromCharCode(65 + i)
                    });
                },
                geocoder: L.Control.Geocoder.nominatim()
            }),
            routeWhileDragging: true,
            routeDragTimeout: 250
        });

var rlayer = L.layerGroup([routing]);
map.hasLayer(rlayer) ? map.removeLayer(rlayer) : map.addLayer(rlayer);

},
    'Display Route'
).addTo(map);

I try use map.hasLayer(rlayer) ? map.removeLayer(rlayer) : map.addLayer(rlayer); but route layer after second click don't hide


texts: Can I compile an addin developed in ArcGIS 10.3.1 for use with 10.1 and 10.2?

texts:
I made an update to one of my addins mainly because users requested new features. I made the changes and compiled the addin in 10.3.1. Although I specified that the addin is compatible with 10.3 and higher I am getting emails from users wondering if I could make it work for 10.2.

ESRI provides a way to save Map Documents for downward compatibility.  Is there a way to compile an addin developed in ArcGIS 10.3 for use with legacy versions (10.1, 10.2)

Note: my ArcObject references are set to non specific version (Specific Version = False)


texts: After upgrade to 10.3.1 Publishing Tools will not start?

texts:
We are running a ArcGIS Server Standard Workgroup on Azure VM.
I just upgraded to ver. 10.3.1 and everything looked fine until I tried to publish anything. See video: http://screencast.com/t/VCtkxfHSY8ZX

It also fails when trying to publish from ArcMap (Error messages: http://screencast.com/t/VBIG6Xv8s0nx)

I installed a fresh Azure VM with latest ArcGIS Server image, and everything looked fine (publishing worked). But as soon as I added Users and Roles to one of the services, the publishing failed (as shown in the video above)

Are there anyone else that has experienced this?
Any ideas on how to fix it?


texts:
Is there a software that can create an animation of a flight/drive path by entering two cities?
I don't want to use the tour function in google earth because I would like to see just the path being drawn from point A to B similar to the animated map function in imovie.
In google earth for example I can't do that and manually recording a tour is very tedious. Any ideas? I need this to be used in a slideshow of pictures, so I want to show first the flight path on a map and then pictures.


texts: Preventing interpolation outside of data extent with LAS Dataset To Raster ArcMap?

texts:
When I use the LAS Dataset To Raster tool in ArcMap I get an unwanted interpolation outside of the actual study area, like this:

According to Esri this can be prevented by changing the definition of the input LAS dataset: http://pro.arcgis.com/en/pro-app/tool-reference/conversion/las-dataset-to-raster.htm

  Including a study area boundary as a clip constraint in the definition of the input LAS dataset is highly recommended. One reason is to prevent interpolation from occurring outside the real data extent of the survey. Secondly, there can be a severe performance penalty when using natural neighbor options if the data area is not properly defined.

But how do you achieve this? I can't find these options in the Las Dataset To Raster menu.

My input values are as follows:

Defining a processing extent doesn't seem to work either.


texts: How to find top 10 greatest values in a GEE image collection?

texts:
I'm trying to categorize temperature extremes for a region. Using the Hycom GEE dataset, I know that using a Reducer I can find the maximum value for each pixel across the image collection (see code below). However, I'd like to find that top 10 or 15 temperatures across each pixel in the collection and average them together to characterize a swath of the hottest days that year. It doesn't look like you can specify how many of the max values you want with ee.Reducer.max(). I'm really new to GEE coding so I'm stumped. How can I get the i number of maximum values?
var dataset = ee.ImageCollection('HYCOM/sea_temp_salinity')
                  .filter(ee.Filter.date('2019-01-01', '2019-12-30'))
                  .filter(ee.Filter.neq('system:index','2017022312'))
                  .filter(ee.Filter.neq('system:index','2019030509'))
var temp_0 = dataset.select('water_temp_0', 'salinity_0');
var max = temp_0.reduce(ee.Reducer.max())


texts:
Can I make a rule in Openlayers2.13 in which I will have multiple properties and values? This is what I have now:

new OpenLayers.Rule({
        name: "NORTH",
        filter: new OpenLayers.Filter.Comparison({
        type: OpenLayers.Filter.Comparison.EQUAL_TO,
        property: "owner", // name of column in table
        value: "NORTH" }),
        symbolizer: {
                pointRadius: "4",
                fillColor: "#FF0000",
                strokeOpacity: ".8",
                strokeWidth: "1"
                }
        })

And I need to have something like this:

new OpenLayers.Rule({
        name: "NORTH",
        filter: new OpenLayers.Filter.Comparison({
        type: OpenLayers.Filter.Comparison.EQUAL_TO,
        property: "column_I","column_II" // name of column in table
        values: ["NORTH","something_else"] }),
        symbolizer: {
                pointRadius: "4",
                fillColor: "#FF0000",
                strokeOpacity: ".8",
                strokeWidth: "1"
                }
        })


texts: Dealing with Null Values in Destination Table after Editing with ArcGIS Relationship Class

texts:
I have an ArcGIS geodatabase that contains one origin feature class linked to multiple standalone tables by Simple one-to-many Relationship Classes. The foreign keys in the destination tables are all nullable, deliberately, in order to allow shapes in the origin feature class to be edited, merged and (in some cases) deleted.  My editing scenario is generally deleting an existing shape and then creating a new one with the same primary key.

ArcGIS behavior sets the foreign key value to "Null" if I delete an origin feature.  Even if I delete one shape and create a new one with the same origin key I have to go back into the many linked tables and re-add the foreign key which got set to Null when I deleted the first shape. I need to be able to delete, split, or otherwise edit the origin features without the foreign key being set to "Null."

Can anyone suggest the best workaround to deal with this problem?

I'm using a file geodatabase on ArcGIS 10.1 and 10.2 (same behavior in each version). One functionality I need to keep is that when I ID the item on an origin feature class, I need to be able to view (and Edit) all the linked tables.


texts: How to tell if the ArcGIS10.3 Image server extension is installed and activated?

texts:
Most extensions within ArcGIS can be activated over the extensions dialog under Tools>extensions. With ArcGIS Server Image Extension I cannot see whether it is installed yet alone activated. I can publish image mosaics directly so it must be installed is some form. How can I confirm this? In  ArcGIS Server Manager under the Site tab the Software Authorization tab
shows no extensions. If I can publish an image mosaic directly as an image service from a registered FGDB, does this mean it has been installed correctly?


texts: How to calculate edge length in Networkx

texts: I want to export a directed weighted graph from a shapefile. I used read_shp function of the Networkx package to export the directed graph which perfectly matches my needs. But I am unable to calculate the length of each edge as line geometries are simplified into start and end coordinates in the output of Networkx. How can I calculate the exact length of the edges (linestrings)?


texts: Raster calculation over multiple layers using a function or do-loop

texts:
I have a bunch of raster calculations that are working, but there must be a more efficient way to do this with a do-loop or function.

Essentially I have a bunch of combinations of: rates (county level, zip level, etc), which correspond with age group, population datasets for each age group, two sets of concentration data sets, and multiple betas and lower and upper confidence intervals for the betas. I am running the below attributable fraction calculation.

y = rate * population * (1-exp(-beta*concentration))

It's not important to know exactly which combinations I have to run, just that there are a lot of them. I have tried if statements that haven't worked, functions using (i in length) which haven't worked, and multiplying stacked rasters (stacks of the rates, stacks of the pollutants, vectors of the betas, etc).

Doing this for every combination, lower and upper confidence intervals, etc, is not efficient and I am bound to make mistakes.

Is there a straightforward way to get started?

INPUTS - all are on the same resolution, extent, cell size.

#Population rasters
adult.pop
all.pop
eld.pop

#Rates
zip.all # Zip, all ages, all-cause, rate per 10,000
cbg.25 # CBG, ages 25-99 years, all-cause, rate per 10,000
cbg.65 # CBG, ages 65-99 years, all-cause, rate per 10,000
co.all # County, all ages, all-cause, rate per 10,000
co.25  # County, ages 25-99 years, all-cause, rate per 10,000
co.65 # County, ages 65-99 years, all-cause, rate per 10,000

#Pollutants
no2.v1  # V1 concentration
no2.v2  # V2 concentration

#Stack rates by age group
rates.adult <- stack(cbg.25, co.25)
rates.all <- stack(zip.all, co.all)
rates.eld <- stack(cbg.65, co.65)

#Betas
beta.adult <- 0.002078254
adult.lower <- 0.000598207
adult.upper <- 0.003536714
elderly.lower <- 0.001105454
elderly.upper <- 0.001365306

CALCULATION

# (1) V1 concentration
elderly.1 <- eld.pop*rates.eld*(10^-4)*(1-exp(-beta.eldery*no2.1))

# (2) V2 concentration
elderly.2 <- eld.pop*rates.eld*(10^-4)*(1-exp(-beta.elderly*no2.2))


texts: What is the appropriate input scale factor for EVI?

texts:
I am trying to calculate the Enhanced Vegetation Index using Landsat surface reflectance data. The formula for EVI using standard constants is:
EVI = 2.5 * ((NIR – Red) / (NIR + 6 * Red – 7.5 * Blue + 1))
Landsat Collection 1, Level 2 surface reflectance data is scaled from 0-10,000 upon delivery, but it is often rescaled from 0-1 for analysis. A simple analysis reveals that EVI is sensitive to the scale factor of the input data, so it is important to scale the input data correctly when calculating EVI.
A simple comparison in R:
#Input data scaled from 0-10,000

library(raster)
landsat_20070621 <- brick("./landsat_20070621.tif")
evi_20070621 <- 2.5 * ((landsat_20070621[[4]] - landsat_20070621[[3]]) / (landsat_20070621[[4]] + 6 * landsat_20070621[[3]] - 7.5 * landsat_20070621[[1]] + 1))

summary(values(evi_20070621))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's
-327.5000    0.5154    0.8160    0.8641    1.1570 1240.0000       385

#Input data scaled from 0-1

landsat_rescale <- landsat_20070621/10000
evi_rescale <- 2.5 * ((landsat_rescale[[4]] - landsat_rescale[[3]]) / (landsat_rescale[[4]] + 6 * landsat_rescale[[3]] - 7.5 * landsat_rescale[[1]] + 1))

summary(values(evi_rescale))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's
-0.1074  0.1414  0.1926  0.2047  0.2526  0.9126     385

Many sources suggest that EVI values are inherently constrained between 0-1, although this is not necessarily the case because it does not have the same normalization structure as a normalized difference index (e.g., NDVI). Nonetheless, Huete et al. (2002) report values in the 0-1 range.
Comparing the values reported in literature to the values above, it seems like maybe an input scale factor of 0-1 is correct, although I cannot find this in writing anywhere. Is there a definitive source indicating the appropriate input scale factor for calculating EVI?


texts:
From a DEM of a watershed, how is possible to calculate the longest flow path?

I know this questions has been asked before but I haven't found a complete solution yet.


texts: How can I use ogr2ogr & sqlite to create a new GeoJSON file that contains only a specific set of rows while renaming two of the keys?

texts:
I need help with a syntax error from a sqlite command in ogr2ogr.
I'm trying to run a query on a geojson file. The query should create a new file that contains only a specific set of rows while renaming two of the keys.
What I have tried so far:
ogr2ogr -f geojson -dialect sqlite -sql "SELECT ID as aid, min as amin, max as amax, FROM contour WHERE amax >4 and amax <10" contour-pt5-middle.geojson contour-pt5-min4-max10.geojson

It results in the following error:
ERROR 1: In ExecuteSQL(): sqlite3_prepare_v2(SELECT ID as aid, min as amin, max as amax, FROM contour WHERE amax >4 and amax <10):
   near "FROM": syntax error

Simply removing a comma before my FROM statement seems to get me a bit closer, but it doesn't carry over the geometry object of coordinates:
ogr2ogr -f geojson -dialect sqlite -sql "SELECT ID as aid, min as amin, max as amax FROM contour WHERE amax >4 and amax <10" contour-pt5-middle.geojson contour-pt5-min4-max10.geojson

creates:
{
"type": "FeatureCollection",
"name": "contour",
"features": [
{ "type": "Feature", "properties": { "aid": 6, "amin": 4.0, "amax": 4.5 }, "geometry": null },
{ "type": "Feature", "properties": { "aid": 7, "amin": 4.5, "amax": 5.0 }, "geometry": null },
...

Here is a screenshot of how my data is structured in contour-pt5-min4-max10.geojson:

Expected output
Should be the same as the screenshot but only including features with "aid" from 6 to 16.


texts: Visualizing irregular temporal data using ArcMap

texts:
I want visualize irregular temporal data which are in structure of below picture:

I have to show the changes of the columns G to J with the changes of time. Some rows have the same spatial coordinates but different time and value dimensions.
I tried to do this using the time slider according to Interpolation layer not moving by time slider using ArcMap, but the interpolation layer could not be displayed temporarily.
I am now looking for a solution that allows me to display data with irregular time intervals in a timely manner.


texts: Using ArcGIS to extend a polygon in one direction up to the edge of another polygon

texts:
I have two shapefiles: the first containing a route defined by a single polygon (black line in the image below), the second containing a series of between 1000 and 4000 small polygons located somewhere inside the larger polygon (green and pink boxes).

I would like to extend each of the small polygons in the direction perpendicular to the route, so that they become long polygons that end at the edge of the outer polygon. These new polygons don't necessarily have to span the entire width of the route polygon, for example where they appear on a bend, as shown below:

I am using ArcGIS 10.4.1 for Desktop with an Advanced license. I've also got access to a Spatial Analyst license. Is what I want possible, and if so how do I do it? I could also use any freely available software other than ArcGIS if that would make this possible.

-----Update-----

Following the suggestion from @Hornbydd and @FelixIP to use Arc's Euclidean Allocation tool (available with a Spatial Analyst license), I now have what I was aiming to produce, as shown in the image below:


texts: Mapbox Studio - View bridges at any zoom

texts:
I'm currently styling a map with Mapbox Studio. My goal is to get bridges over water to be visible at any zoom. I've tried several different approaches so far.

The first and most obvious approach is to work directly with the bridges layers. This works, but the bridges are only visible when the zoom >= 13. There is an option to change this in select data for the layers, but that doesn't seem to do anything.
Using the roads layer. The roads layer also contains bridges within it, but can't distinguish bridges from roads at < 13 zoom. To get the desired result, the roads layer needs to be in front of the water layer, the water layer needs to be in front of the land layer, and the land layer needs to be in front of the roads layer. Land needs to be behind water, because land does not exclude water. This gives me an impossible configuration. No matter which way I try it, something will be out of place.
The only option that I've really found to be somewhat viable thus far is to make the roads layer the same color as the land layer. This then shows bridges over the water and the roads appear hidden on the land.

I'd prefer not to go with my third approach, as I'd like there to be contrast between the bridges and land. I think the most desirable outcome would be using the roads layer, but hiding it while not over water. This way, bridges on land would not appear, but I'm not confident that there's any way to do this.

Is there any way possible to make either of the first two approaches work; or a different approach all together?


texts: Georeferenced raster image not appearing on map in QGIS

texts:
I imported a PNG file into QGIS for georeferencing, and set it to add automatically to map after finishing the process. However, once its job was done and the new layer and joined my list panel on the map, I couldn't see it displayed (all white background still remains).
What could the issue here potentially be? I have checked all its display and opacity settings, so I know the issue isn't there. I opened up the TIF file itself and it was all clear and available there. I've even done this process and gotten through it before, so the georeferencing part on my end shouldn't be the issue. What has gone wrong here, the projection itself, or could the file size be too big? It hasn't loaded for an entire hour already.


texts: Plot ENU points on a map having an origin GPS Lat and Long

texts:
I have a GPS coordinate point to denote my origin:

8.6258,47.3325

Then I have a series of ENU (east, north up) points in meters that represent the distance from my origin:

-10, -1 //East, north (not using up)
-12, -2
-14, -3
-16, -5
....

I would like to plot this on a map, I know with google maps you can simply provide an array of coordinates and it builds it for you:

{[8.6258,47.3325],[8.6258,47.3325],[8.6258,47.3325]...}

I'm not sure how to convert those ENU points to coordinates. Any ideas?


texts: How to add X/Y coordinates and count of overlapping points in a postgresql table?

texts:
How to add to my postgres/postgis table the following columns

X coordinate | Y coordinate | count (number of points in same location)

from an existing table populated with geometries type POINT


texts:
On my Trimble Nomad, I have ArcPad 8 and windows mobile 8. The handheld has been working great up until recently. The machine has been very choppy and slow and when I go to add a point or really do anything I receive the window mobile gwes.exe error.

ESRI help hasn't been able to find a solution. Has anyone seen this before and found a work around or an actual fix?


texts: Migrate in a new programming GIS language

texts:
The last five years I have been working with VBA and ArcObjects for developing a specific application within ArcGIS. That project that began around 2008-2009 carried out with VBA and ArcObjects which is not supported by ESRI anymore. So, I need to migrate in a new programming language and I am not sure which is the best for my case: Python, VB.NET, C# or Java? My aim is to develop specific planning applications within ArcGIS platform. What is the most easiest to learn taking into account that I have experience with VBA and ArcObjects? Which is the most powerful for developing complex GIS applications? Which has the best support?
May I have your opinions please?


texts: Defining bounding box within Add WFS Layer dialog of QGIS?

texts: I use QGIS 1.8 under WIN XP in Standalone Installer Version. I need to limit the extent of the data to be downloaded via WFS because of the large amount of data in the requested layer. In older/different versions of QGIS there used to be a check box within the Add WFS layer dialog to limit the WFS download to the actual view. In my version of QGIS it is missing. What is the best way to set a bounding box for WFS download or how can I get above mentioned check box?


texts:
Is there any way to reproject an individual layer within ArcGlobe? In the "Draped layers" I load the tile raster (file in GeoPackage *.gpkg format) using the "Make raster layer" tool. Then in the "Elevation layers" I add the altitude data SRTM1 (N20E092.hgt) covering the raster.  The problem is that the draped layer in WG WGS_1984_Web_Mercator_Auxiliary_Sphere, and elevation layer in GCS_WGS_1984! Because of this we have a wrong terrain (mountains in the ocean):

The correct terrain is (as an example the OpenStreetMap WMS layer is taken):

Here are spatial references for draped end elevation layers:

For various reasons I can't change the projection in the original GeoPackage data. Is there any way to do this in ArcGlobe? How do I achieve the alignment of layers?


texts: QTiles Plugin maximum store support, intensive computing time?

texts:
We try to compute 187 000 tiles with QGIS plugin QTiles in zoom level 18 (around 3000 square kilometres). There are 25 layers in there (mainly polylines and polygones). It took much longer than usual to compute these (it's not done yet). The first 18 000 tiles took over 12h to compute.

I know all of this depends on computer adjustments, capacity, etc. and it depends on the number of shapes. My concrete questions are:

Is there a maximum memory/store support on the part of the plugin?
Are there other QGIS plugins you can recommend for this task (with these - parameters)?


texts: Selecting features by current map extent in QGIS 3

texts:
I'd like to select all features in the current map extent that also satisfy some other condition via "Select by Expression".
For instance I'd like to do something like this:
intersects($map_extent, $geometry) and  "name" ilike 'Jay%'

Where $map_extent is a geometry that represents the current extent of the map.
I'm using QGIS 3.4.
I've written the following custom expression function which has the desired result but I'm wondering if there is an easier way?
from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args='auto', group='Custom')
def map_extent(feature, parent):
    return QgsGeometry.fromRect(iface.mapCanvas().extent())


texts: How to add a new features to an existing shapefile using Geotools?

texts: I have a shapefile with some points. I would like to add a new points to this shapefile using Geotools.


texts:
I have a simple translation, from GeoTIFF to ECW in FME.  It creates the file just fine, but it does not define the CRS for the ECW.  When loading the ECW into QGIS I have to pick the CRS.  When loading it into ArcMap I get Projection  not defined warning.

In Workbench I have a GeoTIFF reader, a RasterInterpretationCoercer translator to add a transparent alpha channel, and a CsMapReprojector translator.  The source GeoTIFF uses the NAD83(2011) datum, if I reproject to NAD83 the ECW's CRS is defined, but if I set the source and destination to the NAD83(2011) the ECW is undefined.  The CRS is also undefined if I don't have the CsMapReprojector translator on the workbench.

I have tried to set the coordinate system in the ECW writer step but that doesn't work either.


texts: Convert PostGIS geometry + 3 orientation parameters to a KML Camera object using SQL only

texts:
I have some WKB PostGIS Point geometries in WGS84.
I also have some orientation parameters for these geometries, stored as angles in degrees, namely, heading, tilt and roll values.
I want to build a KML for them expressed as a Camera object.
Does it exist a native PostGIS function for that? Or do I have to build my own? I've found that ST_AsKml can do that but for the position only (it seems there is no possibility to also include the orientation):
SELECT ST_AsKML(location) FROM my_point_table LIMIT 1;

which returns:
"<Point><coordinates>-122.4790,37.8110,127.0</coordinates></Point>"

But I do have angular values as well in separate Postgres columns:
heading: 18.0
   tilt: 85
   roll: 0

I wish I could directly output a KML like this using, as shown on the KML webpage linked hereabove, by only Postgres/PostGIS function(s):
<kml xmlns="http://www.opengis.net/kml/2.2">
  <Document>
    <name>Golden Gate tour</name>
    <open>1</open>
    <Camera>
      <longitude>-122.4790</longitude>
      <latitude>37.8110</latitude>
      <altitude>127</altitude>
      <heading>18.0</heading>
      <tilt>85</tilt>
      <altitudeMode>absolute</altitudeMode>
    </Camera>
  </Document>
</kml>

I'm surprised, I must have missed such obvious thing in PostGIS, right?
Or else, I didn't search in the right direction... Hence my question.
Information
These are the only 'ST_AsKML' function's signatures I can see in my schema:

From the Docker image: postgis/postgis:13-3.1
PostgreSQL: PostgreSQL 13.1 (Debian 13.1-1.pgdg100+1) on x86_64-pc-linux-gnu, compiled by gcc (Debian 8.3.0-6) 8.3.0, 64-bit
PostGIS: 3.1 USE_GEOS=1 USE_PROJ=1 USE_STATS=1
By the way, I noticed there are less signatures than in PostGIS 2.4:


texts:
I am working under QGis 3.4.8 Madeira under Windows 10.

I have uploaded two layers on my project.

Layer A is a points layer. Layer B is a polygons layer.

Some points of layer A are located in the area covered by some polygons of the layer B. I want to delete these points.

My final aim is to have a layer with buffer instead of the "free-range" points with the polygons. However, I cannot delete the points overlapping with the polygons.

My question is then:
How to delete SOME points of a layer depending on their overlapping with another layer ?


texts: Difficulties applying zyp.sen() to a RasterBrick object in R

texts:
I have a RasterBrick object that represents a raster time series (say different images of the same scene over a thirty year period) and I'm trying to perform a Theil-Sen regression on this time series.

I'm attempting to do this through the zyp.sen() and calc() functions but I've run into a roadblock.

Here's what I've tried running:

 temp <- list()
 for(i in 1:10) {
     temp[[i]] <- runif(36,-1, 1)
     temp[[i]] <- matrix(temp[[i]], 6, 6)
     temp[[i]] <- raster(temp[[i]])}
 temp <- brick(temp)

 time <- 1:nlayers(temp)

 fun=function(x) { if (is.na(x[1])){ NA } else { m = zyp.sen(x~time)}}

 calc(temp,fun)

However, I keep running into this error:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
cannot use this function

And I'm not really sure why, or what I'm doing wrong.


texts: How to Migrate web maps from ArcGIS Flex to JavaScript API?

texts:
A client has several web maps that were made using ArcGIS Flex. But the client needs all the current web maps to be changed to ArcGIS JavaScript API.

What steps need to be taken to make this happen efficiently?


texts: Oracle VectorLayer without geometry - How to get attributes?

texts:
I have a problem to iterating over a oracle VectorLayer without geometry in QGIS 2.6.

self.uri = QgsDataSourceURI()
self.uri.setConnection("localhost", "port", "xe", "name", "password")
self.uri.setDataSource("", str(tableName), None, "")
oraclelayer = QgsVectorLayer(self.uri.uri(), str(tableName), "oracle")
# check valid
if oraclelayer.isValid():
    QtGui.QMessageBox.information(None, "OK", "The layer is OK!")

attrs = []

# check feature count
oraclelayer.selectAll()
QtGui.QMessageBox.information(None,"featurenumber",str(oraclelayer.selectedFeatureCount()))

# from http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html#iterating-over-vector-layer
iter = oraclelayer.getFeatures()
for feature in iter:
    attrs.append(feature.attributes())
    QtGui.QMessageBox.information(None, "feature", str(attrs[feature]))
QtGui.QMessageBox.information(None, "Datainformation", str(attrs[0]))

The oraclelayer is Valid.
selectedFeatureCount returned the correct number.

After the messagebox "featurenumber" QGIS returned an exception for the messagebox "Datainformation":
list index out of range


texts: Exporting a map layout (both as pdf and image) results in the uppermost y-axis label being partly cut off. It looks like the upper layout boundary is set by the map rather than the label. Any way to fix this?


texts: Measure the distance between points automatically

texts: How can I measure the distance between points referring to the same element in an automated way? In fact I've imported in QGIS buildings downloaded from two different sources and obviously they are not placed exactly in the same position.


texts: Selection of Projection for a small scale project

texts:
I'm quite unsure if this is the right place to ask but I would like to have some advice from experts.

I'm working on a QGIS project in a small area in the south of Portugal. I've been questioning myself about which projection would be the most appropriate for the project.

Since I'm aiming to establish marine protected areas, I believe that choosing one that would preserve areas and keep them as close to reality as much as possible would be the best option (please do correct me if I'm wrong), meaning an equal-area projection.

Can anyone indicate an EPSG that would fit into these parameters?

Until now, I was using the global 4326 and I also experimented with 32629 (WGS 84 / UTM zone 29N)

EDIT: My area of interest has 9 000 km². I would prefer for the error to be as small as possible but I can tolerate if it is not uniform.


texts:
Some GIS maps can be viewed online by using a web browser. We can identify the data-link of the web map by, Right-clicking -> 1) Inspect -> 2) Sources -> 3) Locate the Data Link
Can a web map be imported into QGIS?
Please use the example that is linked below.
Here's a gis map viewer: https://rpgis.isd.lacounty.gov/Html5Viewer/index.html?viewer=GISNET_Public.GIS-NET_Public


texts: QGIS Layer Style - Adding intensity/color if elements overlay

texts:
I have a layer with thousands of points.

I have been asked to represent this layer in way that can be easily seen the level of concentration of points in the layer. So when there is much point concentration, the zone should be darker (or redder i don't mind) and if there are few points the zone should be lighter (or greener).

I don't know exactly how to achieve this but i tried creating a buffer polygon for every point. The problem is that the color in this new layer is not added while elements overlay...

How could i edit the layer's style to achieve so in QGIS?

Another aproach that could be uselful is to obtain a raster image with red zones and green zones depending on the concentraion of points. I prefer to obtain a vector layer but if i can generate both layers better than best.

Regards,


texts: Cache greater than 500mb ArcGIS Server 10.2.1

texts:
Will this affect performance? Why 500? And how can I minimize the size of the cache?
2 service less than 500 MB will be faster than one with more that 500MB?
Thanks


texts: Projecting shapefile of point geometries?

texts:
I am trying to creating a point shapefile that lies within polygon(a feature layer). By creating point geometry objects and filtering the ones that fall within the polygon(feature layer), I created a list (list of point geometry objects as "pointList"). Using CopyFeature_management tool, I created a shapefile out of "pointList". But this doesnot have projection and the attributes has FID, Shape and ID(all 0).

How could I project this point shapefile, the same projection as that of polygon?

Here is what I have so far, (where the extent is calculated for polygon feature layer):

pointList =[]
m=0
     x = random.uniform(extent.XMin, extent.XMax)
     y = random.uniform(extent.YMin, extent.YMax)
     points = arcpy.Point(x, y)
     pointGeo = arcpy.PointGeometry(points)  #creates point geometry object
     polygon = arcpy.Polygon(polygon_array)  #polygon object
     pointWithin = pointGeo.within(polygon)  #gives boolean output.
     if pointWithin == True:
         m = m + 1
         pointList.append(pointGeo)
arcpy.CopyFeatures_management(pointList, outpath)


texts: Draw a grid of polygons in Google Earth Engine

texts:
Let's say that I want to chart daily mean sea surface temperatures per X by X degree cell over time, within some larger patch of the ocean. As per these instructions, this is easily done in GEE using ui.Chart.image.seriesByRegion(). That is, provided you have already defined a FeatureCollection of all the individual polygon regions (i.e. X by X degree cells) first.

Defining this collection of gridded cells is where I am getting stuck. Is there a smart way to automate this step? Perhaps one can even construct the cells from the properties of an imported image collection, e.g. one of the gridded NOAA datasets?

Alternatively, I could draw each of the cell polygons by hand, but this is obviously inefficient and haphazard. I could also write a javascript loop, but loops are supposedly best avoided in GEE as I understand it.

FWIW, in R I would use sp::makegrid() or do something like this.


texts: R circular statistics on a raster time series

texts:
I have a  raster time-series with values from 1 to 36 representing 10 days (dekads) timesteps and NA values.
I need to calculate the average and the standard deviation, taking into account that dekads are circular (36 is last dekad of the year, 1 is the first ones).

Here below is the code for the average (it seems working) but I am not able to write a code for the standard deviation.

Can you check if the code for the average calculation produces the right outputs and help me in writing the ones for the standard deviation.

require(raster)
r <- raster(ncol=50, nrow=50)
s <- stack(lapply(1:length(r), function(x) setValues(r, values = sample(x = c(0:36),size = ncell(r), replace = T))))

# 0 is set to NA

s[s==0]<-NA

#  function
conv <- 2*pi/36  #SET HERE THE time steps

fun1 <- function(m,na.rm=T){
  x1 = Arg(mean(exp(conv*(m-1)*1i),na.rm=na.rm))
  x2 = x1/conv
  x3 = (x2 + 36) %% 36
  return(x3)
}

# MEAN calculation
s_avg <- calc(s,fun1)
s_avg <- round(s_avg,0)
s_avg <- s_avg + 1
s_avg [s_avg==37] <- 1
s_avg

plot(s_avg)


texts: Iterating over layers and exporting them as images with PyQGIS

texts:
In QGIS I have the layer tree with the following layers:

Roads
Transition
Differences
Study_area

I would like to iterate over the layers and export  two different images where Roads and Differences always remain "On". At the end I would like to get these two Images:

Roads and Differences + Transition
Roads and Differences + Study_area

I have tried to use the  script published in the post How to iterate over layers and export them as PNG images with PyQGIS in a standalone script to become the respective Images, the Code that I am using is:

from qgis.core import QgsApplication, QgsMapLayerRegistry, QgsVectorLayer, QgsProject
from qgis.gui import QgsMapCanvas, QgsMapCanvasLayer, QgsLayerTreeMapCanvasBridge
from PyQt4.QtCore import QTimer, QSize
import sys, os

qgisApp = QgsApplication([], True)
qgisApp.prefix="C:\\Program Files\\QGIS 2.18\\apps\\qgis"
qgisApp.setPrefixPath(qgis_prefix, True)
qgisApp.initQgis()

# Required variables with your shapefile paths and names
pngsPath = 'D:/QGIS_environment/'
RoadsLayer = QgsVectorLayer('D:/QGIS_environment/test_datei/Roads.shp', 'Roads', 'ogr')
DifferencesLayer = QgsVectorLayer('D:/QGIS_environment/test_datei/Differences.shp', 'Differences', 'ogr')
otherLayers = {'Transicion': QgsVectorLayer('D:/QGIS_environment/test_datei/Transicion.shp', 'Transicion', 'ogr'),
    'Study_area':QgsVectorLayer('D:/QGIS_environment/test_datei/Study_area', 'Study_area', 'ogr')}
count = 0

canvas = QgsMapCanvas()
canvas.resize(QSize(500, 500)) # You can adjust this values to alter image dimensions
canvas.show()

# Add layers to map canvas taking the order into account
QgsMapLayerRegistry.instance().addMapLayer( RoadsLayer)
QgsMapLayerRegistry.instance().addMapLayers( otherLayers.values() )
QgsMapLayerRegistry.instance().addMapLayer( DifferencesLayer )
layerSet = [QgsMapCanvasLayer(DifferencesLayer)]
layerSet.extend([QgsMapCanvasLayer(l) for l in otherLayers.values() ])
layerSet.append(QgsMapCanvasLayer(RoadsLayer))
canvas.setLayerSet( layerSet )

# Link Layer Tree Root and Canvas
root = QgsProject.instance().layerTreeRoot()
bridge = QgsLayerTreeMapCanvasBridge(root, canvas)

def prepareMap(): # Arrange layers
    for lyr in otherLayers.values(): # make all layers invisible
        root.findLayer( lyr.id()).setVisible(0) # Unchecked
    root.findLayer( otherLayers.values()[count].id()).setVisible(2) # Checked
    canvas.zoomToFullExtent()
    QTimer.singleShot(1000, exportMap) # Wait a second and export the map

def exportMap(): # Save the map as a PNG
    global count # We need this because we'll modify its value
    canvas.saveAsImage( pngsPath + otherLayers.keys()[count] + ".png" )
    print "Map with layer",otherLayers.keys()[count],"exported!"
    if count < len(otherLayers)-1:
        QTimer.singleShot(1000, prepareMap) # Wait a second and prepare next map
    else: # Time to close everything
        qgisApp.exitQgis()
        qgisApp.exit()
    count += 1

prepareMap() # Let's start the fun
qgisApp.exec_()

But I keep getting the message:

  root.findLayer( lyr.id()).setVisible(0) # Unchecked AttributeError:
  'NoneType' object has no attribute 'setVisible'

What could be the problem?


texts: Create a grid in R, change in cell size displaces position of grid

texts:
I'm writing R script which will allow users to overlay a grid on their data in a new projection with the option of defining the size of the grid cell. This seems to work fine, however, when I plot two grids with different cell sizes over top of each other, the grids do not align (see photo). Any suggestions on how to standardize the position of the grid so grids with different cell sizes will align with each other?

new projection to use

projection <- "+proj=aea +lat_1=-30 +lat_2=-50 +lat=-40 +lon_0=175 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"

select size of cell to use

size <- c(250, 1000, 2000, 4000, 8000, 16000, 32000, 64000)

grid.size <- size[6]

read in data

data.df <- "your data set"

turn data into SpatialPoints

data.df <- SpatialPoints(cbind(data$lon, data$lat), proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))

lat/long projection

print(extent(data.df))

plot data in lat/long coordinates

plot(coordinates(data.df))

change lat/long projection to new projection and plot

data.t <- spTransform(data.df, CRS(projection))

data.t.coords <- coordinates(data.t)

plot(data.t.coords, xlim=c((extent(data.t)[1]*1.1), (extent(data.t)[2]*1.1)), ylim=c((extent(data.t)[3]*1.1), (extent(data.t)[4]*0.9)))

plot(data.t.coords, xlim=c((extent(data.t)[1]), (extent(data.t)[2])), ylim=c((extent(data.t)[3]), (extent(data.t)[4])))

define lon0, lon1, lat0, lat1

minx <- min(data.t.coords[,1]) #lon0
maxx <- max(data.t.coords[,1]) #lon1
miny <- min(data.t.coords[,2]) #lat0
maxy <- max(data.t.coords[,2]) #lat1

print(c(minx, maxx, miny, maxy, grid.size))

calculate number of cells in each dimension

Lx <- (maxx - minx)

Ly <- (maxy - miny)

Nbx <- (round((Lx) / grid.size)) #* 1.1 #add 10% buffer to grid

Nby <- (round((Ly) / grid.size)) #* 1.1

build grid using centroid of bottom left cell, cell length and width (equal), number of cells in each dimension; spatialise grid w/new projection

grid <- GridTopology(c(minx, miny), c(grid.size, grid.size), c(Nbx, Nby))

grid <- SpatialGrid(grid, proj4string = CRS(projection))

grid <- as(grid,'SpatialPixels')

grid <- as(grid,'SpatialPolygons')

minx <- grid@bbox[1,1]

maxx <- grid@bbox[1,2]

miny <- grid@bbox[2,1]

maxy <- grid@bbox[2,2]

print(c(minx, maxx, miny, maxy))

plot grid

for(i in 1:length(grid)){
  for(po in 1:length(grid@polygons[[i]]@Polygons)){
    pol <- grid@polygons[[i]]@Polygons[[po]]
    xp <- fortify(pol)[,1]
    yp <- fortify(pol)[,2]
    polygon(xp, yp, col=rgb(0, 0, 0, 0), border='red')
  }
}


texts: Converting 32 bit image to 8 bit palette raster image using FME

texts:
I need to convert 32-bit TIFF image (input) to 8-bit palette TIFF image (output) by using FME.
How can I do this by using palette of output image?


texts: Auto increment features downwards - top of map to bottom map

texts:
is there a method/tool that can auto increment based on spatial location.

I have tried a auto increment by "sort" based on (polygon centroids) X,Y coordinates but this doesn't really work that well if the polygons are irregular in location.

What I want is an unique ID for each property starting from 1 (top of map) to n (bottom of map) but going from left to right and not jumping around.

For example (see below) - I want "14" to be "13", "18" to be "14", "20" to be "15" and so on going down the map.


texts: GeoServer does not show some shapefiles using Add Store

texts:
I am working the GeoServer tutorials from Intro to GeoServer via Packt Publishing.

I'm trying to add a shapefile via AddStore using the GeoServer GUI. I can browse to, "see" and publish shapefiles that come with the GeoServer install and are located in the GeoServer data did, but I can't "see" shapefiles I've downloaded myself, even though I can browse to the directory where they live.

I downloaded a county shapefile from the census site given in the GeoServer beginners guide, the files are consistent with the tutorial, but the directory appears empty when I try to add it via Add Stores.

Do I need to restart GeoServer? Is there a file I'm missing?


texts: Clipping multipatch to polygon using ArcGIS for Desktop?

texts:
I have a multipatch dataset and a 2d polygon.

Which ArcGIS geoprocessing tool can I use to select only those multipatches whoms x,y are inside the polygon?

I want to select all multipatch features that intersect if Z is disregarded.


texts: Convexhull for selected features and performing buffer for the resultant convexhull using python

texts:
I have features selected on the mapcanvas. I want to perform convex hull for those selected features and save those in a memory or store it in a variable preferably and then want to do buffer for the resultant convex hull using python in QGIS.

layer = iface.activeLayer()
features = layer.selectedFeatures()
convex_hull = QgsGeometryAnalyzer().convexHull(layer, path.shp, True, -1, p-None)

The above convexhull function exports the result to a layer but doesnt store. How this can be done?


texts: Rotate feature shapefile with certain pivot point and angle?

texts: I want to rotate a line feature using certain pivot point and angle, but I can not find what tool to use. If rotating a raster, I can simply use Rotate tool, and there I can easily specify the pivot point and rotation angle. But what should I do if I want to rotate a feature? I know I can simply use "Edit Feature" and "Rotate" toolbar, but it does not provide me the way to specify the pivot point and angle.


texts: Advantages either way between GML and Shapefiles?

texts:
After a good start many months ago of getting into GIS through QGIS, I had to drop what I was doing and did not touch it for quite some time. I am now able to get back into it and have a need to load GML vector maps. It strikes me that they have the same properties as Shapefiles, but there must be a difference.

Can someone give a simple explanation of the difference and the advantages either way?

I see that one can be converted into the other, so if there must be a reason to do so.


texts: Problem importing ESRI shapefile information to MySQL

texts:
I recently updated to Qgis 2.0 and therefore the "exporttomysql" plugin does not work anymore. What I want to do is exporting the geometry information of a shapefile(Polygon) into a MySQL database. I found some information on how to that at http://www.bostongis.com/PrinterFriendly.aspx?content_name=ogr_cheatsheet and adopted their code to fit my needs

ogr2ogr -f "MySQL" MYSQL:"CSV_DB,host=myhost,user=user,password=pwd,port=3306" -nln "world" -a_srs "EPSG:4326" ~/Documents/Fields.shp

However, I get the error "Unable to open the datasource ~/Documents/Fields.shp with the following drivers....

When I open the file with "ogrinfo" it tells me the the opening with the driver "ESRI Shapefile" was successful. 1; FieldsFinal(Polygon)

Any ideas how I get the data into my database?

Thanks


texts: ArcGIS Server Tiles as Basemap in CartoDB

texts:
Has anyone successfully used their ArcGIS Server Tilecache as a CartoDB Basemap?  I'm trying to hit mine directly using the XYZ Tile URL Template as follows:

http://mercator10.library.yale.edu/ArcGIS/rest/services/Digital_New_Haven/Atlas_of_New_Haven_1911_Mosaic/MapServer/tile/{z}/{x}/{y}

without success. I'm probably missing something obvious.


texts: Reducing the number of Google Map API calls by using Open Street Map tiles

texts:
I have a website that is exceeding the Google Maps API free limit.  I am wondering if by switching over to Open Street Map as the default tile layer (replacing Google Maps Street layer) if I will reduce the number of API calls?

I am planning on switching over to Leaflet, but that will take longer.

I am paying for each map load. So using Open Street Maps will not help.

  How are map loads counted on the Google Maps Platform?  A single map
  load is charged when any of the following occur:

  A web page or application displays a map using the Maps JavaScript    API.
  An application requests a single map image from the Maps Static    API.

https://developers.google.com/maps/faq#usage_mapload


texts: How to hide Infowindow when using Arcgis Measurement widget?

texts:
I am using infotemplates in ArcGISDynamicMapServiceLayer to show the info popups on the map.

But the problem arises when using Measurement tool the infowindow popups on every click on the map. Since I haven't placed any listeners or code for fetching the features and showing the info window. I don't know how to stop this behavior.


texts: I am trying to create a 1km2 grid of Kenya but I'm struggling. I have changed the projection to EPSG:32337 but this makes Kenya about 10 metres long which is obviously not right and leads to an error that the area I want to grid is too large. The standard projection (WGS84) will only allow me to create a grid in degrees.


texts: Specifying ArcGIS Server WMS GetFeatureInfo tolerance?

texts:
is there a way to specify pixel tolerance when requesting feature info via GetFeatureInfo WMS request for an ArcGIS Server WMS-published map service? I have not been able to find any such information in the documentation.

The OGC WMS standard does not specify any request parameters for coordinate tolerance, I just wonder if there's similar setting somewhere which could make retrieving feature info a bit more convenient for the user, especially for point features.


texts: Why won't Categorized renderer see GeoPackage attributes/values?

texts:
I've encountered a problem with the Categorized renderer in QGIS. I have a GeoPackage file including a set of data points, with approximately 30 attributes (mainly strings). I'm trying to set the colour of points depending on a string attribute value, using Categorized renderer.

When I click on the "Classify" button, the values of the relevant attibute are not loaded and the classification won't succeed: there's only one symbol/colour for all points. The values can't be loaded in the Field Calculator, too: upon clicking the "Unique" or "Sample" button the field beneath stays empty. The Attribute Table on the other hand is displayed correctly, with all values for all attributes in the right cells.

This occurs both in QGIS 2.18 and in QGIS 3.2, in the case of this specific file only. Also, some of the attributes are not affected. The problem can be bypassed by saving the file as Shapefile instead of GeoPackage, however then the attribute labels (Headings in Attribute Table) are cut short, so such approach is not practical.

Could someone help me understand what is going on and how to make it work?

I suppose it may be a text coding issue, but I can't change the coding method for GeoPackage format (the field is greyed out).

The GeoPackage file can be downloaded from here:
https://drive.google.com/file/d/1-GUObpa4yanmX9Y1tBIPAmsh1kgQB9NL/view?usp=sharing

The file contains a list of creches in Olsztyn, Poland, with addresses geocoded using MMQGIS/Google Maps (publicly accessible data only). I've only been able to use the Categorized renderer with the attribute "Telefon" (phone number), but not with other atributes.


texts:
I am a newbie to GIS systems so please be patient.

I have a PostGIS database set up with a spatial data table. One of the columns is of Geometry type and contains Points with SRID 4326. I want to view these points overlaid on a map. I can connect to PostGIS from QGIS, load the points as a PostGIS layer and using the OpenLayers plugin also overlay these on a map. The trouble is, it appears that the latitude and longitudes are reversed when displaying them in QGIS! When I connect to the database from command line and use the SELECT query with ST_ASTEXT on the geometry column, I can see that the lat longs are in correct order.

Any suggestions or pointers on what might be the problem?


texts: How can I get more verbose messages from a GP service?

texts:
When I run a tool I've published as a gp service, the only messages I get are these:

esriJobMessageTypeInformative: Submitted.
esriJobMessageTypeInformative: Executing...
esriJobMessageTypeInformative: Succeeded.

Is there some way I can get the same messages I see when I run the tool locally?

I don't see anything obvious in the documentation.

It is a custom tool written in arcobjects (implementing IGPFunction), so if there's some special way to call IGpMessages I'd consider that.


texts: Downloading surface water monthly multi- images from 2010-2015 in defined region? GEE

texts:
I want to download data of surface water monthly history of region[73,25,105,40]from 2010 to 2015. referencing the others codes I write the codes as follows. But When I Run, there is no action of the GEE.

var AOI=ee.Geometry.Rectangle([73,25,105,40]);
var datasuface=ee.ImageCollection('JRC/GSW1_0/MonthlyHistory').filterDate('2010-10-1', '2015-10-18').filterBounds(AOI);
var list=datasuface.toList(length);
    for (var i=1;i<length;i++){
         var image=ee.Image(list.get(i));
         var name1=ee.String('Monthwater')
                    .cat(ee.String(ee.Number(i)))
                    .getInfo();
         print(name1)
         Export.image.toDrive({
               image: image,
               description:name1,
               scale: 250,
               maxPixels:1e13
         });
    }


texts:
Is there anyway to run an exported script from the processing modeler in PyCharm?

How do I go about installing the necessary packages to run a function?


texts: How to convert a line shapefile into a CSV with nodes and links?

texts:
Is there a way to convert a linear shapefile into a CSV file of the following format:
Node_1, Node_2, attribute_1, ...., attribute_x
     1       2       "hello"        567845.334

I have looked at GDAL/OGR2OGR and PostGIS and I don't think either converts the data to the format I want.
I don't mind writing a script to do it, if I'm guided in the right direction.


texts: Correcting monotonicity errors using ArcGIS Desktop

texts:
There is not a lot of information on this site or any site on how to correct monotonic errors for line vectors in a GIS context. It seems like such a simple problem with a simple fix. Monotonicity is required when a z-enabled feature is required to be ever increasing or decreasing i.e. the z values at each vertex should be consistently going down or up. This is most often seen in hydrology networks than need to be monotonic for network and flow. A monotonic error is flagged when there is a vertex that goes opposite of this trend. See the before and after examples from ESRI below.

Correcting monotonicity errors page from ESRI
Data Reviewer as well as other topology checks can find them. Defense Mapping for ArcGIS claims to fix them by interpolating the z-values before and after the errors to be consistent with the downward or upward slope. After several attempts, The Fix Monotonicity tool in Defense Mapping only makes my lines monotonic by applying one z-value to the entire line essentially flattening it.
Are there any tools or techniques in ArcGIS Desktop that will automatically correct these errors?


texts: Getting text box to stay put in CartoDB?

texts: I'm trying to label two areas of my map in CartoDB. The base map labels two of them already but there are two other regions that need to be labeled. Every time I put a text box over one of the regions, and then move around on the map, the text box moves too.


texts: QuickOSM plugin crashes when using layer extents for query

texts:
I'm trying to use QuickOSM to find some items within the extents of a polygon layer, but every time I run a query with the "extent of a layer" radio box checked, QGIS immediately closes. This happened with every layer I tried, in several different coordinate systems. (Using the map canvas extents worked the way I expected it to.)

I did a little searching, but the closest I came to my problem was this bug report, which was resolved in 2014.

Is there still a problem with layer extents in QuickOSM, or if not, can anyone tell me what I'm doing wrong? Am I using older, or incompatible, versions of QGIS or QuickOSM?

(The versions I'm using are QGIS 2.18.0 on Linux Mint 17.3, and QuickOSM plugin version is 1.4.5.)


texts: Coordinate systems for countries around Caspian Sea?

texts: I'm trying to do spatial analysis of the Caspian Sea area and am not sure which projection to use for the countries around it. The Caspian spans multiple UTM zones and the countries around is are from two continents so I'm hesitant to use one Continental projection. Is it small enough scale picking one of the UTM zones wouldn't distort the countries too much? Just looking for the best projection to use and why.


texts: Why getting 403 error in geoserver 2.10.2

texts:
I know this question asked many times (e.g. OpenLayers 3: Cross-Origin Request Blocked: The Same Origin Policy disallows and
Enabling CORS in geoserver 2.8.2 or how to use JSONP in geoserver?) but those all are not solving my issue.

I am working on WFS_T .Below is my code

    var transactWFS_polygon = function (p, f) {
    switch (p) {
        case 'insert':
            node = formatWFS.writeTransaction([f], null, null, formatGML_polygon);
            break;
        case 'update':
            node = formatWFS.writeTransaction(null, [f], null, formatGML_polygon);
            break;
        case 'delete':
            node = formatWFS.writeTransaction(null, null, [f], formatGML_polygon);
            break;
    }
    s = new XMLSerializer();
    str = s.serializeToString(node);
    alert(str);
    $.ajax('http://xxx:8076/geoserver/xxxx/ows', {
        type: 'POST',
        dataType: 'xml',
        processData: false,
        contentType: 'text/xml',
        data: str
    }).done();
}

I am getting 403 forbidden error in Options.How to enable option in geoserver which is running in localhost.


texts: Automatically fill gaps and merge freshly generated polygons with the longest adjacent polygon

texts:
I get polylines automatically generated from AutoCAD from a colleague. Converting them into polygons with QGIS 3.6, I find out some sliver polygons and gaps in between, as shown in the picture below :

So far, I am correcting mostly by hand which is time consuming, as detailed follow :

To get rid of sliver polygons, I am using the command v.generalize from the GRASS Module, then remove the duplicated geometries and merge manually polygons with the same number.

To fill gaps, I am adding a new feature by clicking on the gap with the tool "Fill gapp with  a new feature" from the digitizing tool plugin, and merge manually polygons.

Any of you have an idea to automate the process, or at least speed it up?

The aim is to get a result like this :


texts: make a selection by attributes of the fourth part of a range of numbers with arcpy

texts:
I am using Arcmap 10.1 and I need to create a script (for use on shapefiles) where I have a input field named "sequence" which is numerical data. I need to select the fourth part of a range of numbers.

For example, 1-100 depending on the amount of data in the shapefile and I need to make a selection by attributes of "sequence" where I divide 100/4 = 25 and the selection is the first  1-25, the second selection of 26-50, the third 50-75, etc.

The number of records can change and be 150/4 = 37 and the selection would be 1-37 the second third etc 38-75

I have the ide of using arcpy.GetCount_management for me to total the records but do not know how to use the resuldato to divide by 4 and then put the result and make the selection with arcpy.SelectLayerByAttribute_management but not how to do this and if this is the right way to do


texts: Replace a field with maximum of it's attribute values based on values repeated on another field's attribute values using ArcPy

texts:
So I have performed Spatial Join between a point and line feature class. Now for each point there are several lines. What I want is the maximum count of line features as an attribute value in the "Count" field based on "label" field same attribute values.
I have thought of few ways to do it but isn't working exactly. I can also add a new field and write the maximum of those "Count" values there.
Here I was trying to convert it into a numpy array and then calculate the maximum.
   '''try:
        arr_label_count= arcpy.da.FeatureClassToNumPyArray("SPATIAL_JOIN_OUTPUT_Splice_Closures_Fibre_Drop_copy",
                                                                     ["label", "Count"], skip_nulls=True)
        arcpy.AddMessage(arr_label_count)
    except Exception as e:
        arcpy.AddMessage("Issue in arr_label_count...: " + str(e))

    try:
        arr_label= arcpy.da.FeatureClassToNumPyArray("SPATIAL_JOIN_OUTPUT_Splice_Closures_Fibre_Drop_copy",
                                                                     ["label"], skip_nulls=True)
        arcpy.AddMessage(arr_label)
        filtered_arr_label = self.remove_duplicates(arr_label)
        # arcpy.AddMessage(filtered_arr_label)
    except Exception as e:
        arcpy.AddMessage("Issue in arr_label_count...: " + str(e))'''

Here is an image of how the array looks -


texts:
I want to reclassify a slope raster as part of Least Cost Path Analysis, but I'm a little bit confused about the reclass rules.

My rules look like this at the moment (numbers left are degree):

0 thru 3.434 = 1

3.435 thru 5.143 = 5

5.144 thru 89 = 20

I'd like to have everything under 3.434 classfied as quite easy to pass through, 3.435 to 5.143 should be harder but still possible and everything over 5.144 basically impossible (or at least very hard to pass). Do the reclassified values represent these thoughts or should I change them?
I also thought about adding more rules to make it more detailed, what do you think?


texts: How to schedule a model to run at a certain time?

texts: I'm trying to schedule a model to run at a certain time and have found this ESRI Blog post which is providing a good start. What I can't figure out though is where in the importToolbox function I tell it which model/tool within the toolbox to run? I have a toolbox that has about a dozen tools/models within it and I just want to run 1 of those.


texts:
I am trying to use ImageMosaic plugin on Geoserver in order to read NetCDF files. I follow this tutorial : https://docs.geoserver.geo-solutions.it/edu/en/multidim/mosaic_config/netcdf_mosaic.html#
When I create ImageMosaic datastore everything is ok.
Here is my indexer.xml :
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<Indexer>
    <domains>
        <domain name="time">
            <attributes>
                <attribut>time</attribut>
            </attributes>
        </domain>
        <domain name="elevation">
            <attributes>
                <attribut>elevation</attribut>
            </attributes>
        </domain>
    </domains>
    <schemas>
        <schema name="ta">
            <attributes>the_geom:Polygon,imageindex:Integer,location:String,elevation:java.lang.Integer,time:java.util.Date</attributes>
        </schema>
    </schemas>
    <coverages>
        <coverage>
            <origName>ta</origName>
            <name>ta</name>
            <schema ref="ta"></schema>
            <domains>
                <domain ref="time" />
                <domain ref="elevation" />
            </domains>
        </coverage>
    </coverages>
    <collectors>
        <collector name="timeregex">
            <value>regex=[0-9]{6},format=yyyyMM</value>
            <spi>TimestampFileNameExtractorSPI</spi>
            <mapped>time</mapped>
        </collector>
        <collector name="elevationregex">
            <value>[0-9]{3}</value>
            <spi>DoubleFileNameExtractorSPI</spi>
            <mapped>elevation</mapped>
        </collector>
    </collectors>
    <parameters>
        <parameter name="AbsolutePath" value="true" />
        <parameter name="AuxiliaryFile" value="_stac.xml" />
        <parameter name="IndexingDirectories" value="/var/lib/geoserver_data/coverages/data/_stac" />
    </parameters>
</Indexer>

My databasestore.properties :
SPI=org.geotools.data.postgis.PostgisNGDataStoreFactory
host=172.17.0.2
port=5432
database=gis
schema=ta
user=docker
passwd=****
Loose\ bbox=true
Estimated\ extends=false
validate\ connections=true
Connection\ timeout=10
preparedStatements=true
create\ database\ params=WITH\ TEMPLATE\=ta

And my NetCDF ancillary file :
 <?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<Indexer>
<schemas>
        <schema name="ta">
            <attributes>the_geom:Polygon,imageindex:Integer,location:String,elevation:java.lang.Integer,time:java.util.Date</attributes>
        </schema>
    </schemas>
    <coverages>
        <coverage>
            <origName>ta</origName>
            <name>ta</name>
            <schema ref="ta"></schema>
            <domains>
                <domain ref="time" />
                <domain ref="elevation" />
            </domains>
        </coverage>
    </coverages>
</Indexer>

And an extract of my PostGIS database table ta :

My problem is that when I want to create a layer my list is empty :

I have no error logs.


texts: Is there a standard color palettes for the spectral indices?

texts: Is there a standard color palette (or most commonly used in publications and other materials) for the spectral indices?


texts: Finding locations within a radius of a given path

texts:
I am trying to figure out in Postgis the right sql statement for the following:

I have in my database a table of merchants with latitudes and longitudes;
I have a route from google that lists a lot of latitude and longitude points for that path;
I need to find all the merchants within 50 miles of that route.

I currently have a query that works for a single point, but it isn't performant to run the query over and over for each point returned by google (there can be thousands). Also, I am using a nodejs server, so I won't be able to multithread the process. Preferably the database would be doing all the work anyways.

Here's the query for the single point:

set search_path to my_schema, public;
SELECT m.*
FROM merchants m
WHERE ST_DWithin(
  Geography(ST_MakePoint(m.longitude, m.latitude)),
  Geography(ST_MakePoint(-94.6855231,38.9287475)),
  ( 50 * 1609.4 ) -- convert miles to meters
);

I don't even know where to start for the other solution.


texts: Openlayers 5: Click event on specific layer?

texts:
I have a couple of layers and i want to add click event on a specific layer of the map.

E.g. If clicking on parcels, do different ajax request than clicking on buildings.

How can i do it?

    var wmsLayer = new ol.layer.Image({
  source: wmsSource
});
var map = new ol.Map({
  target: 'map',
  layers: [baseLayer, wmsLayer],
  view: new ol.View({
    center: ol.proj.fromLonLat([33.33386, 35.14710]),
    zoom: 14
  })
});

map.on('click', function (evt) {
  map.forEachLayerAtPixel(evt.pixel, function (layer) {
    if (layer==='web:Parcels') {
      console.log('parcels clicked');
    }
  });
});


texts: GetPrint Request with SELECTION parameter

texts:
According to informations from this site I can use "SELECTION parameter" with GetPrint Request. I tried to achieve this but without success. With GetMap request everything is OK.

Can some one help me with this?

Here are my tests - URL's and result images.

WMS GetMap request without selection:

http://192.168.1.14/cgi-bin/qgis_mapserv?
map=/home/user/map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetMap&
CRS=EPSG:2180&
BBOX=-1,-1,1,1&
FORMAT=image/png&
LAYERS=test&
WIDTH=500&
HEIGHT=500

WMS GetMap request with selection:

http://192.168.1.14/cgi-bin/qgis_mapserv?map=/home/user/map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetMap&
CRS=EPSG:2180&
BBOX=-1,-1,1,1&
FORMAT=image/png&
LAYERS=test&
WIDTH=500&
HEIGHT=500&
SELECTION=test:0

WMS GetPrint request without selection:

http://192.168.1.14/cgi-bin/qgis_mapserv?map=/home/user/map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetPrint&
CRS=EPSG:2180&
map0:EXTENT=-1,-1,1,1&
FORMAT=image/png&
LAYERS=test&
TEMPLATE=1

And WMS GetPrint request with selection (image look exactly same):

http://192.168.1.14/cgi-bin/qgis_mapserv?map=/home/user/map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetPrint&
CRS=EPSG:2180&
map0:EXTENT=-1,-1,1,1&
FORMAT=image/png&
LAYERS=test&
TEMPLATE=1&
SELECTION=test:0

Here example with FILTER parameter:

http://192.168.1.14/cgi-bin/qgis_mapserv?map=/home/user/map.qgs&
SERVICE=WMS&
VERSION=1.3.0&
REQUEST=GetPrint&
CRS=EPSG:2180&
map0:EXTENT=-1,-1,1,1&
FORMAT=image/png&
LAYERS=test&
TEMPLATE=1&
FILTER=test:%20%22FID%22%20=%200


texts: QGIS: vector layer: changing polygon overlay in an nested polygon

texts:
I have a polygonal layer. The source of the layer is contained in a shapefile.
When I reduce the opacity of the layer, I realize that a polygon "la seine" overlays the polygon "seine et oise"
How to nest them?


texts: How to Use IGraphicsContainer.BringForward?

texts:
I'm making an add-in for ArcMap using ArcObjects.

I'm trying to bring an IElement that I've added to the ArcMap.Editor.Map.ActiveGraphicsLayer layer up one level in Z-Order. How do I add my IElement to a IEnumElement so that I can cast the ActiveGraphicsLayer to an IGraphicsContainer and invoke IGraphicsContainer.BringForward(IEnumElement Elements)?

There is no example documentation for situations other than using the existing user selection.

Update: I've tried doing the VB example in C# like so. The element still doesn't come to the top.

using ESRI.ArcGIS.Carto;

namespace XivicToolbar.ArcMapRelated
{
    internal class MapGraphics
    {
        private static IGraphicsContainer Graphics => (IGraphicsContainer)ArcMap.Editor.Map.ActiveGraphicsLayer;

        public void Indicate(IElement element, int zIndex)
        {
            //Add the element
            Graphics.Reset();
            Graphics.AddElement(element, zIndex);

            //Do z indexing
            for (; zIndex > 0; zIndex--)
            {
                var graphicsSelection = (IGraphicsContainerSelect)Graphics;
                graphicsSelection.UnselectAllElements();
                graphicsSelection.SelectElement(element);
                Graphics.BringForward(graphicsSelection.SelectedElements);
                graphicsSelection.UnselectElement(element);
            }

            //Update the UI
            Refresh();
        }

        private static void Refresh() => ArcMap.Document.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, null, null);
    }
}


texts: Construction Lines offset when added to fabric

texts: We use geo-positioned CAD data then copy paste survey accurate data into our Parcel Fabric.  When we convert them to new parcels, the points move off the correct location by two to 5 feet.  What causes this offset and how can we keep this from happening?


texts: Why Editing a GeoPackage table with ogr is very slow?

texts:
I want to edit some fields in a GeoPackage (download a.gpkg)

# -*- coding: utf-8 -*-
from osgeo import ogr
from time import clock

print ("Start",clock())
source = ogr.Open("c:/.../a.gpkg",update=True)
layer = source.GetLayerByName( "263c6845-e2c3-4eee-ae26-f0e9539f3c2bP" )
laydef = layer.GetLayerDefn()

i=0
feature = layer.GetNextFeature()
while feature:
    i=i+1
    # make some stuff
    # .....
    layer.SetFeature(feature)
    feature = layer.GetNextFeature()
source.Destroy()
print ("Finish",i,clock())

  ('Start', 3.4136520136791866e-07) ('Finish', 1612, 243.56666043106233)

It works, but it takes a long time (4 minutes for 1.600 records).

What is going wrong?


texts: clipping stacked polygons with onderlying polygon layer

texts: I'm currently performing an network analysis where I calculated catchment areas for health facilities in a city. These catchment areas are overlapping and I want to clip these catchment areas with an underlying population layer which is devided in multiple polygons to see the individual amount of people living in each catchment area. If i use the 'clip' tool in qGIS, it merges all these catchment areas but that is not what I want(see image). I want to keep the individual catchment areas with their unique ID, clipped with the underlying layer without merging results. I currently doing it manually(by deleting all the catchment areas except for one), but that takes me ages. Is there a faster way?


texts: Export ArcGrid Raster to Esri file geodatabase on MacOS

texts:
I have downloaded this dataset, what I believe is an ArcGrid Raster, from this site. However, I ideally would like the data to be in a file geodatabase format so that I can access it using GeoPandas in Python. I am on a Mac and don't have access to a Windows machine, so I don't have the ability to export this in ArcMap as the tutorial that I found here suggests.
Is there a way I can convert this without using ArcMap / using software on MacOS?


texts: Adding multiple layers with timed delay in Mapbox/Leaflet

texts:
So I've got a Mapbox map showing data split into several layers. I'd like the user to be able to click a play button to call the function "gogogo", and then the layers get added one at a time every 1 second or so.

I'm trying to use simple Javascript to iterate through an array of addLayer's to add them one at a time, but the result of this code is that all the layers are added at once when the user clicks the button. Is there another way I should do this? Or is there some way through Mapbox or Leaflet's JS libraries to natively slow down the addition of layers?

var map = L.map('map').setView([57.8915,86.8359], 3);

// this just adds the basemap only
addLayer(L.mapbox.tileLayer('mymapboxaccount.7xu07ldi'), 'Base Map', -1);
function addLayer(layer, name, zIndex) {
    layer
        .setZIndex(zIndex)
        .addTo(map);
}

var animation;

function gogogo() {
  animation = setInterval(addNextLayer, 1000);
}

var layerAdder = function(layer) {
    layer
        .addTo(map);
};

function addNextLayer() {
  for (var i = 0; i < arrayLength; i++) {
    layerAdder(layers[i]);
  };

  var layers =  [addLayer(L.mapbox.tileLayer('mymapboxaccount.wv29be29')),
      addLayer(L.mapbox.tileLayer('mymapboxaccount.d08p8pvi')),
      addLayer(L.mapbox.tileLayer('mymapboxaccount.0tivygb9')),
      addLayer(L.mapbox.tileLayer('mymapboxaccount.lv8l4n29')),
      addLayer(L.mapbox.tileLayer('mymapboxaccount.mlr9hpvi'))
      ];

  var arrayLength = layers.length;

  clearInterval(animation);
}


texts: What is the most efficient way to display 30 features in Openlayers and having them update every 1 sec?

texts:
I have data in PostGIS of GPS coords and I want to show in OL 2.13 in real time (the GPS device transmit data to the server every 1 second) - also have Geoserver if needed.

My current solution is this and I was wondering if there is better way ?

 var coords = new OpenLayers.Layer.Vector("coords", {
    projection: epsg4326,
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
        url: http://xxx/yyy/getGeoJson.php,
        format: new OpenLayers.Format.GeoJSON()
    }),
    styleMap: new OpenLayers.StyleMap({ 'default': {
        externalGraphic: 'resources/images/${iconName}',
        graphicWidth: 175,
        graphicHeight: 173,
        graphicOpacity: 0,
        label: '${gpsName}',
        fontColor: 'black',
        fontSize: '12px',
        fontFamily: 'Courier New, monospace',
        fontWeight: 'bold'
    }
    })
});

and I just redraw the layer each second.

Is there better way / correct way to do such a thing ?


texts: Resolution of NED (National Elevation Dataset)

texts:
I downloaded some DEMs from the USGS with a resolution of 1/3 degrees for different areas in the US.

Now I am looking at the files and saw that all DEMs have the same number of columns and rows.

I measured the dimensions of the area covered by different DEMs and found that all have a similar length in N-S direction, but that their dimension in W-E direction varies a lot, being larger towards the south.

So now I am a bit confused. I know that the length of one degree longitude is larger near the equator than at the poles, so if my DEMs have a resolution of 1/3 degree I thought the number of columns in the DEMs for areas located in the north would be smaller than for areas located towards the south, but this is not the case.

On the other side, if DEMs at the north and south have the same number of columns, but cover areas of different dimensions, then the resolution must be different in both cases ...

Can someone help?


texts: Combining Layers not outputting correctly in QGIS?

texts:
I have two shapefiles I'm trying merge together. I want to use this shapefile in another application I use. When I import them in QGIS, they look exactly how I want them to be positioned (relative to one another). The yellow dots in the image below are from one file (Utah city locations) and the purple dots are from another file (Wyoming city locations)

The CRS for the Utah cities is "EPSG:26912, NAD83 / UTM zone 12N" and the Wyoming ones is a generated one (which I assume is because there isn't a projection file for that shapefile: "+proj=lcc +lat_1=41 +lat_2=45 +lat_0=41 +lon_0=-107.5 +x_0=500000 +y_0=200000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

The problem is, when I attempt to merge them using the vector data management tools (taskbar: Vector-->Data Management Tools-->Merge Shapefiles to One), I get a result that looks like this:

I have no idea why this is happening, I would think I should get the output that I'm seeing visually before I merge. Here's a screenshot of the merge tool:

Any ideas on how to fix this issue?


texts: Method in QGIS to do a select by radius but ignore "islands"

texts:
I am looking for a process in QGIS to run a selection from a given point based of a radius that ignores "islands". For example you may have a body of water which acts as a natural barrier past which I do not want to select any features.

Please see my image as an example:

As you can see in this example the selection was made from a point on mainland England but polygons on the Isle of Wight has been selected even though they are separated by the Solent.

It may also be the case that the multi polygon itself does not contain an island but the selection creates one (see above).

I am not necessarily looking for a single tool or plugin because I don't think it exists, but is there anything I can incorporate into the graphical modeler to get around this? I'm thinking of a rule that ignores any polygons that do not touch.


texts:
after a while I have updated my QGIS to 2.18.15 and I have noticed that some saga tools run very slowly. Some months ago I used the saga multilevel b-spline interpolation tool and it runs very fast. But now It needs more than two times to get the same result. Did something change regarding saga? I have the feeling that QGIS needs time to start SAGA in the background.

Any thoughts?


texts: Correct use of SelectLayerByAttribute_management in ArcPy?

texts:
What are displayed in the tableName are polygon features in a geodatabase. I would like for each row in a loop to take each row TableName and apply the sql query from the sqlstring field and make the feature according to that query that corresponds.

    qname   TheTableName SqlString
0   shp_q1  AST_A        SELECT ...
1   shp_q2  AST_B        SELECT ...
2   shp_q3  AST_PX       SELECT ...

I want to make a new feature polygon according to what the query describes, named as the qname column is named in that row.

Example: Take the AST_A feature polygon, apply the query that it has in the SqlString column and make a new feature polygon based on that query called shp_q1

Here I access the columns that are needed:

datasetList = arcpy.ListTables("*")
for dataset in datasetList:
    with arcpy.da.SearchCursor(dataset, ["TheTableName","SqlString"]) as cur:
        for row in cur:
            query = row[1] #to access the query part of the row
                arcpy.SelectLayerByAttribute_management(row[0], "NEW_SELECTION", query)

row[0] is the 'TheTableName' column and row[1] is the SqlString column.

How should I change the SelectLayerByAttribute_management function parameters to achieve that?


texts: GeoServer GetFeatureInfo - how to get only one result per layer

texts:
I am using GetFeatureInfo to create a popup for different layers and features in a map in Leaflet. I have up to 3 layers that can be visible and overlapping at any given time. If I put the option "feature_count=3" in my request to GeoServer it works great with the exception of when I am zoomed out and features are very small. When this happens and only one layer is visible I will often get 3 results in the popup for that layer. It seems that since the location I click on might refer to more than one feature in a layer, it returns three features for that layer. I have different code that inserts only visible layers into the code below.

How do I ensure that when clicking the map, only one feature's info per layer is returned in the popup? I wish the "feature_count" option referred to how many features to return per layer.

'http://xx.xx.xx.xx/geoserver/workspace/wms?SERVICE=WMS&VERSION=1.1.1&REQUEST=GetFeatureInfo&LAYERS=layer1,layer2,layer3&QUERY_LAYERS=layer1,layer2,layer3&STYLES=&BBOX='+BBOX+'&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&feature_count=3&X='+X+'&Y='+Y+''

Is there an option in GeoServer that I am missing or an option in my code? Do I have to generate a different link like the one above for each layer?


texts: Detecting duplicate attribute values in table using ArcPy without Advanced/ArcInfo level license?

texts:
In ArcMap, the right-click Summarize option very quickly creates a frequency listing of the chosen field:

Note that the count field shows the count of duplicate values. How can I access this using ArcPy?

Is it possible to access ArcGIS's table operations "Summarize..." and "Statistics..." in python? suggests using Summary Statistics, but I can't get any meaningful output from this tool. Instead, it always just lists the total count of features, rather than the frequency of each individual value:

How can I achieve the result from the second screenshot above, using ArcPy in ArcMap 10.0 with an ArcEditor license (whereby I can't run the FREQUENCY command)?

Dan Patterson's script applies to ArcGIS 9.3 so (probably) isn't useful?


texts: Using rasters with different dimensions in GDAL_calc.py

texts:
I am processing several large (>40GB) geotiffs with GDAL_CALC.py and I encountered this error when attempting the process:

  Error! Dimensions of file image.tif (198000, 126000) are different
  from other files (198752, 126000). Cannot proceed.

This was the command I use to replace problem areas in a DEM with values from another DEM:

  gdal_calc -A problem_areas.tif -B mosaic_A.tif -C fill_mosaic_B.tif
  --outfile=dem_corrected.tif --calc="((A!=999)*B)+((A>=999)*C)"

My rather convoluted solution is to produce a convex hull vector from one raster and then use that to clip it and all other rasters with gdalwarp, thus redefining the dimensions:

  gdalwarp -q -cutline convex_hull.shp -tr 0.5 0.5 -of GTiff
  -srcnodata 0 -dstnodata 0 mosaicA.tif mosaicA_clipped.tif

I imagine each clip will take at least a day to process (one is processing slowly as I write this). Does anyone know if it is possible to calculate (using GDAL_calc.py) rasters that have different dimensions without resampling or clipping the rasters? Is there an easier way to avoid adjusting the dimensions, or is this the correct (albeit painful) approach? For instance, can the one problematic file be clipped to the dimensions of the other DEMs?


texts: Invalid SRID in PostGIS when importing data by ogr2ogr

texts:
I am using PostgreSQL 12 with PostGIS 3.0.0 and GDAL 3.0.0 to copy data in it. The command I am using is:

ogr2ogr -f "PostgreSQL" PG:"host=127.0.0.1 dbname=db user=user
 port=5432 password=1" "file.geojson" -nln mutable

The issue here is when i check the SRID of the output table it is 900914 which is not correct in my case it should be  2223. To fix this i have to manually specify EPSG code by adding '-t_srs EPSG:2223'. This was not needed with previous version of ogr2ogr and PostGIS 2.5.


texts: Wrong dimensons on .png-file when using gdal_translate

texts:
I am using

gdal_translate -of PNG -ot UInt16 -scale MIN MAX 0 65535 "source.tiff" "output.png"

in the command-prompt to create a .png-file from a .tif-file that I created using the Interpolation-plugin. But when i view the png-file in Windows photo, it has the wrong dimensions, it gets really distorted(squished) in the y-axis, and less noticeably distorted in the x-axis. But if I drag the .png to QGIS it isn't distorted anymore.

This is the output of gdalinfo for the .tif:

Size is 3000, 3000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,298.257222101],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (-142.699000000000012,251.647000000000048)
Pixel Size = (0.068380000000000,-0.119760000000000)
Corner Coordinates:
Upper Left  (    -142.699,     251.647) ( 10d30'35.92"E,  0d 0' 8.17"N)
Lower Left  (    -142.699,    -107.633) ( 10d30'35.92"E,  0d 0' 3.49"S)
Upper Right (      62.441,     251.647) ( 10d30'42.54"E,  0d 0' 8.17"N)
Lower Right (      62.441,    -107.633) ( 10d30'42.54"E,  0d 0' 3.49"S)
Center      ( -40.1290000,  72.0070000) ( 10d30'39.23"E,  0d 0' 2.34"N)
Band 1 Block=3000x1 Type=Float32, ColorInterp=Undefined
  NoData Value=-9999

And for the .png:

Size is 3000, 3000
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,298.257222101],
            EXTENSION["PROJ4_GRIDS","@null"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",15],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
Origin = (-142.699000000000012,251.647000000000048)
Pixel Size = (0.068380000000000,-0.119760000000000)
Corner Coordinates:
Upper Left  (    -142.699,     251.647) ( 10d30'35.92"E,  0d 0' 8.17"N)
Lower Left  (    -142.699,    -107.633) ( 10d30'35.92"E,  0d 0' 3.49"S)
Upper Right (      62.441,     251.647) ( 10d30'42.54"E,  0d 0' 8.17"N)
Lower Right (      62.441,    -107.633) ( 10d30'42.54"E,  0d 0' 3.49"S)
Center      ( -40.1290000,  72.0070000) ( 10d30'39.23"E,  0d 0' 2.34"N)
Band 1 Block=3000x1 Type=UInt16, ColorInterp=Gray
  NoData Value=-9999