texts: When a class in an ArcGIS extension should be an ArcGIS class

texts:
Visual Studio templates provided by the ArcGIS SDK offer a "ArcGIS Class" which is a simple C# class that has the attributes Guid, ClassInterface, and ProgId before the class declaration.
While adding a new class to an ArcGIS extension project, when should I use an ArcGIS class and when should I use a regular C# class?


texts: Seeking Australian Irrigation Boundaries?

texts: Does anyone know of a source for Irrigation district boundaries in Australia, specifically NSW?


texts: Grouping ZIP code boundaries to delineate territories in GIS

texts:
I am using ArcGIS 10.5.

I have a list of ZIP codes and a number of listings (consider it the number of incidents in each polygon). In each ZIP code I want to group the ZIP codes in such a way that I get territories (bigger polygons by merging the shape) having a total of listing (sum of incidents in all the polygons merge for making the bigger polygon) in range 700-750 per territory.

I tried the grouping analysis tool, but I could not apply the logic of bringing listing value to 700-750 in each territory.

Steps:

Take a polygon if the listing value (attribute) is less than 750. For that polygon look for another nearest polygon to it and merge it with the first polygon.
Add the listing attribute of both polygons to get the value for the bigger polygon.
Keep on repeating 1-2 till the bigger polygon reaches the 700-750 range.
Continue for the whole map.


texts: Adding field in one shapefile which lists where each attribute of that file overlaps with another shapefile?

texts: I have two shapefiles, one has points where people divert water from rivers, the other has polygons of land ownership. I would like to add a field in the attribute table of the points shapefile which lists which land ownership polygons each point overlaps with.


texts: Carto time-series widget: how to change bucket labels?

texts:
I am using the Carto builder animated time-series widget based on a number column that contains a year as a four digit integer.  In the widget's horizontal axis, the buckets are not labeled with integer values but with decimal values (e.g. 1.8k instead of 1800).  People using the map are confused by this.
How can I configure the time-series to use four digit integers for the labels in the timeline?


texts: Changing column datatype with QGIS Field Calculator

texts:
I have seen a solution that I don't get to work - I need to change data type from string to integer and I did what is shown i the image below
Can someone point out what I do wrong?


texts: Wrong length values on QGIS 3.0.3 and 3.2

texts:
While using QGIS 3.x today I noticed that length values were not calculated and measured correctly on EPSG:32611. Calculating and measuring the same data on QGIS 2.18 latest LTR and 2.2 on EPSG:32611 does provide the correct length values.

I already tried on two different computers and the results are the same (wrong values on 3.x, correct values on 2.x). Also importing data generated on other sources (AutoCAD and Draftsight) and calculating length values provides the same result.

I know there are similar questions already asked but none have satisfactory answers. Is anyone else having similar issues and can provide help?

Steps to reproduce:

Create new project on QGIS 3.x on EPSG:32611 (or maybe other UTMs).
Create new shapefile line layer on EPSG:32611 (or maybe other UTMs).
Create a line object and use advanced digitizing tools to force a length value (5).
Use field calculator and/or measuring tool to get the length value of the newly created object.
Repeat steps 1 through 4 on QGIS 2.x using the same EPSG code and compare length values.
(optional) Import data from other sources on QGIS 3.x and 2.x, calculate/measure length values and compare.

EDIT: per Aaron request here is one of the similar questions already asked Field Calculator returning incorrect length and area values


texts: Add elevation point data and lines to topography plot using ggplot2

texts:
Hopefully this is a simple query. I have pieced together some code to make a topography plot using elevation data and ggplot2.

However, now I would like to overlay onto this plot the country outlines, along with points to show the location / elevation of points of interest e.g. cities and also other polygons and lines to show for example routes. However I am struggling to see how best to do this. Any help or suggestions would be appreciated.

My code is as follows:

library(dplyr)
library(ggplot2)
library(raster)
library(scales)
library(rgeos)

#get elevation data
dem1<- getData("SRTM",lat=34.5,lon=33)
dem2<- getData("SRTM",lat=34.5,lon=36)
dem <- merge(dem1,dem2)

#create polygon to crop the elevation data file
library(sp)
coords = matrix(c(34, 31,
                  34, 32.5,
                  36, 32.5,
                  36, 31.0,
                  34, 31),
                ncol = 2, byrow = TRUE)
P1 = Polygon(coords)
Ps1 = SpatialPolygons(list(Polygons(list(P1), ID = "a")), proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

#crop the elevation data using the polygon
dem <- crop(dem,Ps1)

#lower the reolution to enable faster plotting
dem_lower_res <- aggregate(dem, fact=10)
dem.p  <-  rasterToPoints(dem_lower_res )
df <-  data.frame(dem.p)
colnames(df) = c("lon", "lat", "alt")

#plot the data
ggplot(df, aes(lon,lat)) +
  geom_raster(aes(fill = alt))+
  scale_fill_gradientn(colours = terrain.colors(10))


texts: Python Error: TypeError: Append() takes at most 5 arguments (6 given)

texts:
I am writing a python script and have set my variables. But i keep getting this error
"TypeError: Append() takes at most 5 arguments (6 given)"

I had a look online and it says restart my python window but that doesnt work.
Any ideas. Below is a simple version of what i have.  The append code works in Arcpro so i am assuming my setup is incorrect in IDLE.
from datetime import datetime
d = datetime.today().strftime('%b%y') #'Dec20'

# Local variables:
i.e
Leases_x = L:\\{}\\SA\\production\\product\\leases.gdb\\Leases".format(d)
Leases_Y = L:\\{}\\SA\\production\\product\\leases.gdb\\LeasesY".format(d)
Leases_Z = L:\\{}\\SA\\production\\product\\leases.gdb\\LeasesZ".format(d)
... ect

# Process: Append
arcpy.management.Append("'Leases_x ';'Leases_Y '; 'Leases_Z'; 'Leases_S'; 'Leases_T'", 'Leases_A', 'NO_TEST', "", "", "")


texts:
apt-get install python-gdal installs version 1.6.0 and this is already installed. I would like to install the 1.9.0 version.

Do I have to de-install 1.6.0 first and then compile manually or is there an apt-get package that I have missed?


texts:
MSSQL has a great feature called Table Valued Parameters. It allows you to pass a table of a custom data to stored procedures and functions.

I was wondering what is the equivalent , if exists, in PostgreSQL, using JDBC? I know about the option of passing Arrays as a function parameters, but that seems limited to Postgresql data types.

Consider the following PL/pgSQL code:

CREATE  TYPE number_with_time AS(
_num   float,
_date  timestamp
);

and this function header:

CREATE OR REPLACE FUNCTION myfunc(arr number_with_time[])

Can anyone post a JAVA code using JDBC driver of calling that function with an array of the user defined data type?

Thanks


texts:
First of all I would like to say that I am new to gis-se and as such I will try to comply as much as I can to your policy. That said I am actually doing a little project which as the title says relates to RTK-GPS. I am trying to make use of RTKLib which would be familiar to the people reading this.

Following the paper written by the creators of RTKLib I am replicating their low cost RTK receiver with of course more up-to-date parts such as U-blox LEA-6T instead of the U-blox LEA-4T. However, as I am in the UK I do not have access to a free RTK base station within 20 Km range. I therefore looked for a way of making a base station which from a logical deduction could be a second receiver (fix) acting as a base station (for double differenciating) and connected to the Rover (first moving receiver). Anyway, after a brief research I found this setup that seems to be good enough.

My questions so far are the following:

Would that setup work with my receiver?
Is it good enough? Can I still achieve cm-level precision?

If you any other source of information or examples that would help get a better grasp of RTK-GPS I would be very appreciative.

During my research, I read about the whole theory and did understand it to an acceptable level. However, I still have one thing that I cannot understand. For RTK, we perform an ambiguity resolution which consists in finding the integer cycles of the wavelenght. There are many techniques (even some working with floats) but why an integer number gives a better precision? Is it because of the PLL's? Because the wavelength for L1 for example is about 19 cm so I guess an integer number would give a precision within 19 cm, is that correct?


texts: "Length" field needs "ft" added to end of each number!

texts: I have a field for "Length" regarding some data collected in the field. Unfortunately the "ft" after the number for the length had spaces or no spaces or a weird format, but basically I want Length to look like this:             "9 ft" -     "7 ft"  -   "8 ft"  ...etc. At the moment it just says, "9, 7, 8." Anybody know how to fix this in the field calculator in ArcGIS Pro?


texts: Georeferencing CIA's Lambert Conformal Conic projection in QGIS?

texts:
I'm trying to load into the QGIS an image generated from the NA CIA Fact Book PDF  which uses Lambert Conformal Conic Projection, standard parallels 25°N and 77°N. I'd like to work with the exact projection from this PDF and overlay some vector data on top of it.

I think there is something wrong in my workflow. Here is what I've done so far

export the PDF to a TIFF file
Define a Custom CRS in QGIS:

  +proj=lcc +lat_1=25 +lat_2=77 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs

Open the tiff file in georeferencer.  When prompted for CRS, pick the custom defined lcc as above
Add 6 or so GCP at the grid intersections
Transformation settings: Linear/Nearest Neighbor/Target SRS: select the custom defined as above.
Once the image loads into qgis it does not look correct. I've tried Polynomial 2 transform but that didn't fix either.

My understanding is that when the source and dest projection does not change I should be able to use the Liniar transformation.

using QGIS 2.18.

result in qgis:


texts: GeoServer WPS process fail with OpenLayers 2: "The request body should be contained in a CDATA section"

texts:
I am trying to create an application with WPS builder like this: http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/wps.html
When input format = "text/xml; subtype=wfs-collection/1.0", it creates input field for name of the layer on my server:

function addWFSCollectionInput(input) {
var name = input.identifier;
var field = document.createElement("input");
field.title = input["abstract"];
field.value = name + " (layer on the server)";
addValueHandlers(field, function() {
    input.reference = field.value ? {
        mimeType: "text/xml; subtype=wfs-collection/1.0",
        href: "http://46.28.111.137:8080/geoserver2711/wfs",
        method: "POST",
        body: {
            wfs: {
                version: "1.0.0",
                outputFormat: "GML2",
                featureType: field.value
            }
        }
    } : undefined;
});
document.getElementById("input").appendChild(field);

}

This input is required for example by algorithms gs:Clip and gs:Aggregate. When I'm trying to run some of these algorithms, it returns this kind of result:

<?xml version="1.0" encoding="UTF-8"?><wps:ExecuteResponse xmlns:xs="http://www.w3.org/2001/XMLSchema" xmlns:wps="http://www.opengis.net/wps/1.0.0" xmlns:ows="http://www.opengis.net/ows/1.1" xmlns:xlink="http://www.w3.org/1999/xlink" xml:lang="en" service="WPS" serviceInstance="http://46.28.111.137:8080/geoserver2711/ows?" version="1.0.0"><wps:Process wps:processVersion="1.0.0"><ows:Identifier>gs:Clip</ows:Identifier><ows:Title>Clip</ows:Title><ows:Abstract>Clips (crops) features to a given geometry</ows:Abstract></wps:Process><wps:Status creationTime="2015-08-25T08:35:59.568Z"><wps:ProcessFailed><ows:ExceptionReport version="1.1.0"><ows:Exception exceptionCode="NoApplicableCode"><ows:ExceptionText>The request body should be contained in a CDATA section, otherwise it will get parsed as XML instead of being preserved as is</ows:ExceptionText></ows:Exception></ows:ExceptionReport></wps:ProcessFailed></wps:Status></wps:ExecuteResponse>

How can I fix it?


texts: make a loop that can select segments in polyline one by one

texts: The goal is to write a python script that will allow me to move on a polyline making a selection segment by segment as shown in the pictures.  I hope to find a loop that can make that !!


texts: Filter Sentinel-2 images base on cloud cover over a region of interest

texts:
I know that there are a lot of examples of how to mask clouds for sentiel-2 images, but what I want is sort the images over an area of interest without change anything of the image.

Applying a cloud mask in zones where there are to much clouds, changes dramatically the image (in this case in over Bogota-Colombia). And, in some cases it generates holes in the image. thats why i dont want to change anything about the image.

I am taking this code Filter Landsat images base on cloud cover over a region of interest
as reference.

and here is how it will be with Sentinel-2 image.

//the parameter for sentiel-2 images for clouds is
//CLOUDY_PIXEL_PERCENTAGE
var ic =ee.ImageCollection('COPERNICUS/S2');

// A polygon representing the roi.
// my_polygon es a random polygon over Bogota Colombia (very cloudy zone)
var geometry = my_poligon;

var c = ic.filterBounds(geometry);

var withCloudiness = c.map(function(image) {
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
  var cloudiness = cloud.reduceRegion({
    reducer: 'mean',
    geometry: geometry,
    scale: 30,
  });
  return image.set(cloudiness);
});

var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
print(filteredCollection);

//my image with less clouds over my zone of interest
var my_image = ee.Image(filteredCollection.filterBounds(centroid)
    .filterDate('2019-01-01', '2019-06-30')
    .sort('CLOUDY_PIXEL_PERCENTAGE')
    .first());


texts: Copy features from feature class to an existing feature class with arcpy

texts:
I see CopyFeatures_management, but this creates a new feature class.

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000035000000

Is there a way to copy features from a feature class into an existing feature class?


texts: Seeking file representing extents of coordinates systems and their wkid?

texts:
I'm looking for a file that represents the major coordinate systems throughout the world, and the extent that they cover.

Ideally they'll also have the wkid as an attribute.

Something like the shapefile found here of UTM zones, but with wkid's
http://earth-info.nga.mil/GandG/coordsys/grids/universal_grid_system.html


texts: Generate coordinates with Minimum/Maximum distance from given coordinates

texts:
Background:

I am trying to create a coordinates generator for Unit Testing my code.
I would like to try to test cases for methods like:

// Check if within 10km.
// Should pass because it is 0 - 10km from location
location.isWithin(generateLocation(location, 10), 10) === true

// Check if within 10km.
// Should fail because it is 20 - 25km from location
location.isWithin(generateLocation(location, 25, 20), 10) === false

I have tried lots of Stackoverflow answers I never found one answer to this problem.

What I have tried:

function generateNearbyLocation(latitude, longitude, max, min = 0) {
  // Convert max in kilometers to degrees
  const radiusInDegrees = (max * 1000) / 111300;

  const u = Math.random();
  const v = Math.random();

  // random radius
  const w = radiusInDegrees * Math.sqrt(u);

  // random angle
  const t = 2 * Math.PI * v;

  const x = w * Math.cos(t);
  const x1 = x / Math.cos(latitude);
  const y1 = w * Math.sin(t);

  // Not working attempt to add Minimum distance
  // const newLatitude = latitude + (min * 1000 / 111300) + y1;
  // const newLongitude = longitude + (min * 1000 / 111300) + x1;

  const newLatitude = latitude + y1;
  const newLongitude = longitude + x1;
  const distance = getDistanceFromLatLonInKm(latitude, longitude, newLatitude, newLongitude);

  return {
    newLatitude,
    newLongitude,
    distance
  };
}

This works but is not accurate. Even with max of 10, it sometimes generate 11km coordinates.

Here's a playground that I made: https://jsbin.com/hakemubogo/edit?js,console

EDIT

I found a much more accurate formula from another answer but still can't create minimum distance.

function generateLocation(latitude, longitude, max, min = 0) {
  const EARTH_RADIUS = 6371; // km
  const DEGREE = EARTH_RADIUS * 2 * Math.PI / 360 * 1000; // 1° latitude in meters

  const r = (max * 1000) * Math.random() ** 0.5;
  const theta = Math.random() * 2 * Math.PI; // random * (360deg in radians)

  const dy = r * Math.sin(theta);
  const dx = r * Math.cos(theta);

  const newLatitude = latitude + dy / DEGREE;
  const newLongitude = longitude + dx / (DEGREE * Math.cos(deg2rad(latitude)));
  const distance = getDistanceFromLatLonInKm(latitude, longitude, newLatitude, newLongitude);

  return {
    newLatitude,
    newLongitude,
    distance
  };
}


texts: shortest distance between road and lake with PostGIS

texts:
I have a country-wide dataset with lakes (polygons) and streets (lines). For each lake I want to know the shortest distance to the closest street. (I do not necessarily need the closest street itself ... )

I use this query for it:

SELECT
  *,
  (SELECT
      ST_Distance(lake.geom, street.geom)
    FROM
      street
    ORDER BY
      lake.geom <-> street.geom
    LIMIT
      1
  ) AS distance_closest_street
FROM
  lake

This query works quite fast on regional level, but takes forever on country level.  I have the usual "gist"-index on the geometry columns.

Do you have ideas how this query could become faster?

Edit:
I use PostgreSQL 9.5.7 and PostGIS 2.3.2


texts: Enable layer plugins when using external standalone QGIS python scripts

texts:
I created a QGIS standalone script (custom Application) that enables me to automatically export Atlas as images.

The problem is that in my QGIS project, I'm using a Mask plugin and the generated images don't take into consideration this plugin (even if I install it on /usr/share/qgis/python/plugins) when executing my script in a standalone way whereas when I manually export Atlas as images using QGIS desktop and the same project, the plugin is activated and the mask is done correctly in the generated images.

I also tried to use the API sys.path.append('/usr/share/qgis/python/plugins') just after the qgs.initQgis() and also QgsApplication.setPrefixPath("/usr/lib/qgis", True) (tested also with /usr/share/qgis andn /usr folder) before the QgsApplication([], True) without more success.

Is there a specific API to use when executing standalone scripts so that QGIS engine enables layers plugins?

I also created an issue on the Mask plugin project in case.


texts: Google Earth Engine, imagecollection aggregate_max (or min, or min), return null. What's the problem?

texts:
I am a little confused about this function in ee.ImageCollection:

  aggregate_max(property)

  Aggregates over a given property of the objects in a collection, calculating the maximum of the values of the selected property.

It returns NULL, and I do not know why. Here is the code:

// creat a function for clipping.
var clippedtoGulf = function(image) {
  return image.clip(geometry);
};

var ndvi = ee.ImageCollection('MODIS/006/MOD13Q1')
  .map(clippedtoGulf)
  .filterDate('2010-01-01', '2017-03-01');

print('ndvi', ndvi)
Map.addLayer(ndvi, {}, 'ndvi')

var stats = ndvi.aggregate_stats('NDVI')
var max   = ndvi.aggregate_max('NDVI')
print('stats', stats)
print('max', max)

https://code.earthengine.google.com/2ecef30aa9c1294af95eaaa53dccdf06


texts: My team is working on an application that has a map feature, and we want to extend this feature greatly. Among other things, we want the user to be able to load a CAD drawing, and place it on the map. The map is just a single image right now, but we want a more advanced and flexible solution, like eg. tile based. Some of our users have showed us some CAD drawings (DWG) they want to use, and they are absolutely huge, so we are afraid that performance will be bad. I have snooped around and found that (eg.) bing maps support map tile overlays, and I guess that converting the huge CAD drawings to overlayed map tiles could make performance much better, while still maintaing a certain level of control, like clicking layers on and off. But how can I convert the CAD drawings to map tiles from within our application? Is there any API out there that can help with that? In C#?


texts: Calculating area fraction of pixels covered by larger pixel using Google Earth Engine

texts:
My goal is to calculate an area fraction from one image as the threshold to select (mask)another image.

Specifically, the land cover layer has a spatial resolution of 30m and the EVI data is 250m. I want to calculate the area fraction of corn coverage (using the land cover layer) for each EVI pixel. Then the EVI pixels whose fraction values are larger than 0.7 will be extracted.

pixel-based area fraction for corn = (area of corn) / (area of each EVI pixel)

Although the GEE has the reduceRegions method, it requires to input features rather than rasters, and I did not find a method to do this task.

Could anyone help me?

Here is the data:

// import evi and crop data layer (CDL)
var evi = ee.ImageCollection("MODIS/006/MOD13A2")
            .filterDate('2015-01-01', '2015-12-31')
            .select('EVI')
var cdl = ee.Image("USDA/NASS/CDL/2015");

// mask other types of crop and only keep the corn
var corn = cdl.select('cropland').eq(5).updateMask(cdl.select('cropland').eq(5))


texts: Raster Calculator Python Error 00539 for friction surface

texts:
I am trying to convert a slope raster to a friction surface with values based on Tobler's Hiking Function in ArcGIS Pro, but I get an error message I am not sure is about the formula syntax or file/path names. I first converted a DEM of part of the Florida coast with a 5-meter cell size and elevation units in centimeters to a slope raster with a Z-factor of 0.01 because I believe I had to convert centimeters (vertical units) to meters (horizontal units). When the Z-factor was 1 most of the land, including the coastline, had values greater than 45 degrees. Other than roads and reservoirs the slope raster is very flat but I would expect that for Florida.
Error:
 ERROR 000539: Traceback (most recent call last):
  File "<expression>", line 1, in <module>
  File "<string>", line 8, in rcexec_DC684390_0C54_4F14_9EBF_B34DF3F2E1B1
RuntimeError:  ERROR 010240: Could not save raster dataset to cm5_fl7_hike.tif with output format FGDBR.
 Failed to execute (RasterCalculator).

Parameters:
Map algebra expression:  (5.0 / 1000) / (6 * Exp(-3.5 * Abs(Tan(("DEMs\5M\cm5_fl7_slopem.tif" * 3.14159) / 180) + .05)))

Output raster: d:\FLGIS\Lidar\Lidar.gdb\cm5_fl7_hike.tif

I am trying to replicate the following mathematical formula which was typed exactly as it was written in an article guide


texts: Overlay Analysis in QGIS (or Preferably Python)

texts:
I would like to perform an overlay analysis in Python or QGIS to output a list of each shapefile that falls within certain boundaries. I have unsuccessfully tried to use the Spatial Join, Intersect, and Union commands that were described in this similar GIS Exchange post.

For each parcel described in this shapefile (links to download), I would like to see whether it exists within the regions of this shapefile.

Would you please point me in the right direction on how to do this? If there is a Python library that I can use, I would prefer that method.


texts:
I am using the ArcGIS API for Javascript and am wondering if it comes with functionality for zooming out using a rectangle, similar to the behavior in ArcGIS Desktop?  Specifically, the user draws a rectangle and it zooms out based on that rectangle.

I can create the functionality myself, but was wondering if it already existed in the API...


texts: Writing For Loop with Group of Rasters for Batch Conditional Test?

texts:
How do I ask a “for loop” to operate on groups of three files folder (e.x. files A, B, C), then files B, C, D, and then files C, D, E, and so on?

I have a folder full of rasters (.tif): RasterA.tif, RasterB.tif, RasterC.tif, RasterD.tif…and so on. I’ve got a python script that can run the batch conditional test I want. However, it only runs through one raster at a time, and I’d like it to run through multiple.

I think it’s easier to show via code. Please see below.

# Omitted the top part of the script for readability.
# I have all the imports (arcpy and os), set the environment, set overwrite, set cell size
# Please let me know if I need to add this top part back in to the question and I will

rasters = arcpy.ListRasters("*", "ALL")
for raster in rasters:
    inRas = Raster(raster) #pulls the first raster (ex. Raster A)
    inRas2 = Raster(raster)
    # Right now I know this means inRas and inRas2 are the same thing.
    # What I'd like to do is have inRas2 be the next raster "in line" (ex. Raster B)
    inRas3 = Raster(raster)
    # Same deal as above, would like it to be Raster C

    # The CON test itself
    output = arcpy.sa.Con(inRas >= 25, 1, 0) + arcpy.sa.Con(inRas2 >= 20, 1, 0) + arcpy.sa.Con(inRas3 >= 15, 1, 0)

    # My output raster needs to be identified by the first file in the group
    # This part runs fine, can ignore in context of question
    RasName, RasExt = os.path.splitext(raster)
    output.save(outWorkspace + RasName + ".tif")

The fixed script is below. Everything ran without problem.

# Get a list of Rasters
rasters = arcpy.ListRasters("*", "ALL")
# set it so can pull group of rasters with indexing
n = len(rasters)
# run the con test on batch, in groups of 3
for i in range(2,n):
   rasterA = Raster(rasters[i-2])
   rasterB = Raster(rasters[i-1])
   rasterC = Raster(rasters[i])
   output = arcpy.sa.Con(rasterA >= 25, 1, 0) + arcpy.sa.Con(rasterB >= 20, 1, 0) + arcpy.sa.Con(rasterC >= 15, 1, 0)
   RasName, RasExt = os.path.splitext(rasters[i-2])
   output.save(outWorkspace + RasName + ".tif")
   # I have it print just to make sure naming files goes right
   print RasName
   print RasExt
   print "Done"


texts: How to create circles starting with points and attributes (diameter, lat/long)?

texts:
I have a point shapefile with attributes (diameter, latitude  and longitude). For each point in that file, I want to create a polygon (circle) described by those attributes. The new shapefile should have circles of described diameter.

I am using ArcGIS 10.0.


texts: Are there any intelligent travelling salesmen?

texts:
Jokes aside, I have had a routing problem that is almost a travelling salesman problem (TSP):

the starting point is defined
the end point coincides with the starting point
each node has to be visited
the total cost should be minimised

Two years ago I thought TSP would be a perfect match, so I ran some sample data through tsp_solve and Concorde. Luckily, it was quickly obvious that the TSP shortest path is not the real shortest path, since the problem is made easier by unrealistically requiring the nodes be visited exactly once. This picture is just a one-step manual attempt at optimisation of the computed solution and already it saves about the distance of the longest used edge.

The problem resurfaced, as I'm trying to find optimal routes to subsets of mapping/monitoring sites. Location and road network data is both pretty accurate and precise, so an exercise like this makes sense.

I've looked at generalisations of the TSP, but didn't find an appropriate algorithm. Minimum spanning trees don't account for returning from branches (the 1st solution here costs 3 more). From what I understand, the shortest path problem eventually only cares about two nodes and those out of the optimal path would be left out. A special case of the vehicle routing problem seems to fit best, though I don't know if it considers non-direct paths.

My question: is there any settled name, definition for this kind of problem (family)? What algorithm and tool would you use to solve it?

I'm sure it would be computationally heavy, but I'm interested in both general (infinite resources) and practical answers.


texts:
I started using the GDAL Library in a c++ project. It is a great library, but I am a bit stuck right now with different aspects of it.
What is the best way to use gdal_contour in my c++ programm? Should i use the system() function or simply grab the source and use it?
I tried to get the contours of a hightmap with gdal_contour as a PNG or SVG(or different vector format, like eps or so) but had no success so far. Is there a list wich output formats are supported. Until now I only got shp and dxf files.
Do I have to burn the shp file to a png with gdal_translate? I tried this and had some troubles with the alpha channels. Maybe someone has the commands for me to get a transparent contour as png. Sometimes the contour lines look a bit to pixelated, is there a way to smoothen things up with more antialiasing or so?
Is there a way to get a colored contour automatically, or do i have to use a color table txt file?

I tried the following where output.tif was a transparent tif. but had no success.

gdal_contour -i 50 <DEM> contour
gdal_rasterize -burn 0 -l contour contour output.tif


texts: How to call Legend end points of multiple layers published from Geoserver?

texts:
I have map services published from ArcGIS server being called from Angular- Leaflet app and we are displaying the legend (controlled through map controls) from
https://pxxxx.com/arcgis/rest/services/aaa/aaa/MapServer/Legend?f=pjson&token=ddddd.&callback=angular.callbacks._0

and receiving the Json as below based on layerID
 angular.callbacks._0({
 "layers": [
  {
   "layerId": 1,
   "layerName": "abc",
   "layerType": "Feature Layer",
   "minScale": 0,
   "maxScale": 0,
   "legend": [
    {
     "label": "",
     "url": "76d5be60e5691543860aa46358f5dcf8",
     "imageData": "iVBORw0KGgoAAAANSUhEUgAAABQAAAAUCAYAAACNiR0NAAAAAXNSR0IB2cksfwAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAF9JREFUOI1jYaAyYIExzFcw/EeXPBnBwEisPIaBlACYZScjGBjhBmKzjZBrsAGquBDDQGQnU8VAagKqGIjss8HpQgwDKY0MrOmQWmA0HZIH6J8OKSkbqZIOCbqQkrIRAPKcIWF3tuu1AAAAAElFTkSuQmCC",
     "contentType": "image/png",
     "height": 20,
     "width": 20
    }
   ]
  },
  {
   "layerId": 3,
   "layerName": "xyz",
   "layerType": "Feature Layer",
   "minScale": 750000,
   "maxScale": 0,
   "legend": [
    {
     "label": "",
     "url": "0c13bfd1080ed2b63ddc93e033c6739e",
     "imageData": "iVBORw0KGgoAAAANSUhEUgAAABQAAAAUCAYAAACNiR0NAAAAAXNSR0IB2cksfwAAAAlwSFlzAAAOxAAADsQBlSsOGwAAAElJREFUOI3t1DEKACAMA8AIfi+vzQPrVEFwqpmkmTpdM2XCnAkAkuIVIjk26IikIDkOML9UsLxtDRtssMGPQccu2hoeA1vdwVsWVeQROSDODxUAAAAASUVORK5CYII=",
     "contentType": "image/png",
     "height": 20,
     "width": 20
    }
   ]
  },

The most complex map service consists of 53 layers and there are some more services. We are planning to switch to geoserver. Geoserver's getlegendgrahics works only for one layer at one time. So almost feels like there has to be 53 + calls just to gather legend info (unless I am missing something). Also legend in current application are being displayed based on the max and min scale through map controls- which in case of geoserver are stored in SLDs- is there any workflow similar in geoserver e.g. what will be the end point to use SLD info? (it may sound as 2 questions but both are from same context and very much related.)


texts: How to import OpenStreetMap data into GraphX?

texts:
How to import OpenStreetMap data into a Spark GraphX RDD?

Would like to see feasibility of doing drive-time analysis in our own Spark cluster.

Thought to use osm2pgsql or a similar project, but not sure if preserves data required for routing.

Looked at OSM2Routing, but it's not clear what config file required e.g. for whole U.S. map?

I'm new to both OSM and GraphX.

Update. GraphX needs two tables. 1st - Vertex Table (id of a node, properties); 2nd - Edge Table (SrcId, DstID, properties). https://spark.apache.org/docs/latest/graphx-programming-guide.html#example-property-graph


texts: Modifying OpenLayers 2 to connect to different WMS using JavaScript?

texts: I have a program that runs in tomcat and uses JavaScript to access a postgreSQL database as well as Google, Yahoo, Geoserver, and Open Street Map map services through OpenLayers 2.12. I would like to modify the parameters so it will connect to a different WMS; I would like some guidance on where to start looking for these parameters, woulod they be in the OpenLayers.js?


texts:
I'm looking to do spatial joins of billions of points to get the city name for each point globally using spark. I'm not able to find a good publicly available shapefile/geojson/kml of city polygons to pull this from.

My current thought is to source this from open street maps. I've already extracted administrative boundaries and all place=city,town using osmosis and imported it into PostGIS. It looks like place=city,town is not populated consistently and there are major cities missing with this tag, so I'm not sure its usable (e.g. [48.860959, 2.352662] is in Paris, but there is no place=city).

So, I'm looking at resorting to using administrative boundaries. The trouble I'm having with administrative boundaries is that different countries define cities at different admin_levels. Even within the US, NYC is admin_level 5 and every other city in the US seems to be 8. Is there a simpler way to pull this or do I have to create a lookup table with an entry for each country that determines what levels the city comes from in preference order? For example:

US=5,8
GB=6,8
DE=6,8
FR=8
...

If so, does this lookup table of country to admin_level for cities exist anywhere?


texts:
Is there an easy way to get .gpx files pulled from my Garmin GPS60 with MapSource into MapInfo 12?

I know I can grab the data as .txt and fiddle with it in Excel first, just wondering if there is another way.

Thanks.


texts: R gdalwarp issues under new PROJ6 libraries

texts:
i have a workflow that incorporates gdalwarp from gdalUtils in R, to reproject 1km UK data in British National Grid to 0.01 degree data in WGS84. This is always a tricky transformation due to the awkwardly similar cell sizes, but i found that gdalwarp, using the cubic transformation, struck the best balance in terms of maintaining surface total and also cell values;
# i know that sp and CRS is deprecated here but has had little effect on projectRaster etc.
BNG <- CRS("+init=epsg:27700")  # British National Grid
LL <- CRS("+init=epsg:4326") # For LatLong

# BNG file being read in to transform
class      : RasterLayer
dimensions : 1350, 980, 1323000  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : -230000, 750000, -50000, 1300000  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs
source     : C:/Path/BNG1km_file.tif
names      : my_BNG_raster
values     : 0.0001111782, 873.1203  (min, max)

# structure of WGS84 raster being trasnformed to is in the gdalwarp function but for clarity;
class      : RasterLayer
dimensions : 1250, 1840, 2300000  (nrow, ncol, ncell)
resolution : 0.01, 0.01  (x, y)
extent     : -13.8, 4.6, 49, 61.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : memory
names      : layer
values     : NA, NA  (min, max)

# gdalwarp function
gdalwarp(paste0("C:/Path/BNG1km_file.tif"), dstfile = paste0("C:/Path/LL0.01_file.tif"),
s_srs = BNG, t_srs = LL, r = "cubic", tr = c(0.01, 0.01),
te=c(-13.8, 49, 4.6, 61.5), output_Raster=TRUE, overwrite=TRUE, verbose=TRUE)

Checking gdal_installation...
GDAL version 2.1.1
GDAL command being used: "C:\Program Files\QGIS 2.18\bin\gdalwarp.exe" -overwrite  -te -13.8 49 4.6 61.5 -tr 0.01 0.01 -s_srs "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs" -t_srs "+proj=longlat +datum=WGS84 +no_defs" -r "cubic" -of "GTiff" "C:/Path/BNG1km_file.tif" "C:/Path/LL0.01_file.tif"
class      : RasterBrick
dimensions : 1250, 1840, 2300000, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.01, 0.01  (x, y)
extent     : -13.8, 4.6, 49, 61.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : C:/Path/BNG1km_file.tif
names      : layer

# new surface
class      : RasterLayer
dimensions : 1250, 1840, 2300000  (nrow, ncol, ncell)
resolution : 0.01, 0.01  (x, y)
extent     : -13.8, 4.6, 49, 61.5  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
source     : memory
names      : layer
values     : -1621.413, 7279.047  (min, max)

Last year this exact function gave me a good reprojection in terms of data  totals and values, but this year it went haywire (21 times the surface total). Also when using the function against last year's data again, it also produced very poor results.
Has gdalUtils not taken into account the GDAL and PROJ6 updates? or has gdalwarp gone wrong somewhere?
edit: i can see the gdalwarp help (GDAL documentation, not R) has a -ct argument for WKT strings in the update, but this has not made it to R yet
https://gdal.org/programs/gdalwarp.html


texts: Heatmap still renders to concentration rather than attribute

texts:
I have house price data for the UK plotted by postcode, rendered as a heatmap.

I have styled the heatmap to weight points by transfer price and not concentration of points, however I'm running into issues. (See below)

Highlighted point contains 12 houses (with same postcode) at average house sale of £196,721

Highlighted point contains 8 houses (with same postcode) at average house sale of £192,375

Highlighted point contains 3 houses (with same postcode) at average house sale of £194,666

So with all three areas containing the same average house price, the clusters are still weighted by concentration of points.

However, there are times when the heatmap seems to work:

The settings:

What the heck is going on?


texts: How to change position of columns in attribute table?

texts: If I have an attribute table with columns A B C D, is it possible to move the columns to different positions, say, B C A D? When making a new column with the field calculator it would be useful to put it somewhere other than last position or to be able to move key columns close to it.


texts: How to change fieldName value in identify popup panel

texts:
I'm newbie in this ArcGIS API for javascript.
The problem is I want to change the fieldName value . In SQL Query, it is like

Select ....,Length_m_,  CASE WHEN Length_m_ <= 0 OR
Length_m_ IS NULL THEN '-' ELSE CONVERT(nvarchar(50), CONVERT(decimal(10, 8), Length_m_)) END AS newLength.

The output is like this

Here is my code:

/// Bridge 5: {
title: 'Bridge',
description: [
'
',
'
.........	',
'',
'
',
..........
'
',
'Length(km):	',
'{Length_m_}	',
'
',
'
',
'Lane:	',
'{Lane}	',
'
',
'

'
].join(''),
}

So the output will be like;

How I can solve this porblem?


texts: Selecting/highlighting point feature in ArcGIS Server using JavaScript

texts:
How do you go about "selecting/hightlighting a feature" in the map control when a user clicks on a point in the map, all via ADF JavaScript?

I know it can be done server side with a selection graphics layer etc, so I am just wondering if it can be done client side.


texts: Find adjacent cells on raster and compare altitude

texts:
I would like to obtain a list of neighbouring cells for raster data. Below is a simple and intuitive way of doing it using the adjacent() from the gdistance library. This function however only allows finding neighbouring cells for one specified cell. To obtain neighbouring cells for all, I use a nested loop which is horribly inefficient for large data sets. How could I do this more efficiently? I thought about apply but have not managed. The goal is to compare the altitude of each cell on the grid.

library(raster)
library(gdistance) # package for adjacent function

# make the blank raster, the extents and resolution of which are up to you.
# the cells are numbered 1 to ncell (starting in the top left in R)
r <- raster(xmn=0,xmx=6,ymn=0,ymx=6,res=c(2,2))
r[] <- 1:ncell(r)

adj_cells <- c() # initialise vector to store adjacent cells
adj_cells_temp <- c() # initialise second vector to store adjacent cells for each row on the grid
system.time(
  for(i in 1:r@nrows){
    for(j in 1:r@ncols){
      adj_cells_temp <- as.data.frame(adjacent(r, cells=c(i, j), directions=8, pairs=T))
      adj_cells <- rbind(adj_cells, adj_cells_temp)
    }
  }
)

The result also includes duplicates which is not helpful to increase efficiency.


texts: Calculating length of lines using PostGIS?

texts:
I have a set of lines in a table, I want to calculate the length of each line in kilometers and add is a new column in the table. How would I do this?

my table has been set up as

CREATE TABLE public.line (id integer Not Null,
                          name VARCHAR (25),
                          geom geometry (LineString, 27700));

when I try to calculate the distance of the line in km I get an error message saying only lon/lat coordinate systems are supported in geography. how do I fix this?


texts: Item Description on a Feature Class Not Appearing on a Published Feature Access Service via ArcGIS Server?

texts:
I have a feature class in an SDE where metadata has been loaded via the 'Load XML Recordset Document' option of the feature class. When viewing the item description
of the feature class, the metadata appears (which is what I want), such as:

After I publish the feature class as a Feature Access Service via ArcGIS Server, then view the item description of the feature class from the service. The metadata does not appear, such as:

Any thoughts to have item description on the feature class appear on the Feature Access service?

I am using ArcMap 10.3.1, ArcGIS Server 10.3.1, and a SQL Server SDE.


texts: Is there a way to get centroid of geographic object in MSSQL?

texts:
I have a table of OSM geographic roads in MSSQL 2012.

I want to create a new column with the centroid of each road (So I can set center on OpenLayers map for each road).

How can I do it ?


texts: Restrict text widget to only integer between 0 and 999 in IntraMaps Roam

texts:
I have created a form and added an attribute as a Text widget in the Roam manager.

Now I'm trying to restrict the input to be an integer between 0 and 999. I want the field to stay "empty" (or "blanc") as long as the user enter a letter and only to accept integers (but not more that 3 digits). If a 4th digit is entered, the field should continue only showing the 3 first allowed digits.

In QGIS I was able to achieve this with QIntValidator (but QIntValidator is not "available" i Roam?):

global gangeField
gangeField = myDialog.findChild(QLineEdit,"GANGE")
# gangeField will only accept inntegers between 0 and 999
validator = QIntValidator(0, 999, myDialog.findChild(QWidget,"tab2enhet"))
gangeField.setValidator(validator)
gangeField.textChanged.connect(felt_onTextChanged)

I have tried to add some Python (ver. 3.3.2) code in the __init__.py file, but only with limited success:

    def accept(self):
    """
    Called before the form is accepted.  Override this and return False to cancel the accepting of the form.
    :return:
    """
    # if ( int(self.boundwidgets['gange'].value()) == False ):
        # self.pushError("Angi integer!")
    if ( int(self.boundwidgets['gange'].value()) > 999 ):
        self.cancelload("Angi tall mellom 0 og 999")
    return True

Alternative: Can I change the error message to view my personalized text (e.g."Gange must be a number between 0 and 999.") for 2 seconds (like the Warning that is visible for 2 seconds if a required field is missing)?


texts: Elevation Profiles Using QGIS vs. Google Earth

texts:
I have several routes I am calculating elevation stats for.  What I want is Total Ascent/Cumulative Gain, where all ascents and descents are accounted for.  As opposed to solely the difference between the min and max elevations. https://en.wikipedia.org/wiki/Cumulative_elevation_gain

I've been using the QGIS plugin Climber, which provides an output like such

Where Climb (total ascent?) is just maxelev-minelev + descent.  I've been able to verify these numbers by exporting the z values to spreadsheet and using the same formula.

My dilemma is when I check these paths in Google Earth: the profiles, max, and min are all the same, but the elevation gain/loss appear wildly different. For example, for the fourth record above Google Earth provides a gain/loss of 369/99.

Is there a something conceptually I am not understanding?

Edit:
Both profiles:excuse the reversed directions and feet/meter unit discrepancy.


texts: segmentAlongLine() returning line with a gap in the middle

texts:
I'm trying to use segmentAlongLine() to return a portion of a polyline based on percentage in my feature class. It works for most of my lines in the feature class, but in some cases, the polyline returned contains a gap in the middle.
For example, here is the original polyline.

When I create a new line using polyline.segmentAlongLine(0, 0.5, True), I get the following line as an output:

Has anyone else run into a similar issue or know what would be causing the gap? It happens for any value of end_measure greater than 0.27 (27%). I've also tried switching the start_measure and end_measure values but end up with the same issue.
I'm using Python 2.7 with ArcMap 10.7.1.
EDIT: one detail I failed to mention is the polyline is the result of a Dissolve with "DISSOLVE_LINES" set as my unsplit_lines attribute.


texts: Does anyone know if I can use ArcMap or ArcCatalog to change the projection of a .sid img without changing it to a .tiff? I have tried the project raster tool but it changes it to a .tiff which will not be usable for the outcome I am looking for.


texts: Small coordinate offset on QGIS using CRS ETRS89 / UTM zone 30N

texts:
I'm currently trying to produce simplified field slips on QGIS for an undergraduate mapping project in Villamanin, Northwest Spain. I've imported vector map layers from the download centre from the IGN for Spain and have been working through them to produce slips suitable for mapping contacts/etc onto. When starting to render slips I've noticed I have a small coordinate offset between raster maps downloaded from Iberpix and the corresponding slip I've produced on QGIS (see image). Does anyone have any ideas about this?

Whilst the offset is small I'm concerned it will introduced real errors into my final project map...


texts: How to make QGIS show raster image with correct axis scale? NOT SQUARE!

texts:
I've loaded two .IMG raster files from USGS into QGIS. The IMG files themselves are perfect squares of about 10,800x10,800 pixels at 10meter resolution.

But in reality, the terrain they represent is NOT square because the pixel scale is 8 meters for xScale and 10 meters for yScale.

This is confusing because QGIS is displaying the images such that they look perfectly square. How can I make QGIS scale the image in each axis independently so that I get a true visual representation of the map/terrain proportions?


texts: How to fix rasterio error: The EPSG code is unknown. Unable to open EPSG support file gcs.csv

texts:
I try to open a jp2 image with the Python's Rasterio package. The code is below. The problem is that the coordinate reference system is not recognised and when I want to see the meta attribute of the opened band I see that the crs is None.

import rasterio
b6_path = 'D:\\IKI\\1np\\myimage_band06.jp2'
b6 = rasterio.open(b6_path, driver='JP2OpenJPEG')
b6.meta

#output:
# {'driver': 'JP2OpenJPEG',
# 'dtype': 'uint16',
# 'nodata': None,
# 'width': 5490,
# 'height': 5490,
# 'count': 1,
# 'crs': None,
# 'transform': Affine(20.0, 0.0, 699960.0,
#        0.0, -20.0, 4900020.0)}

As I know what the coordinate system is, I tried to set it manually like this:

import rasterio.crs
crs = rasterio.crs.CRS({"init": "epsg:32634"})
b6.crs = crs

But when I call the CRS() function I get this error message:

CRSError: The EPSG code is unknown. Unable to open EPSG support file gcs.csv.  Try setting the GDAL_DATA environment variable to point to the directory containing EPSG csv files.

I checked the GDAL_DATA variable by this command in the Anaconda prompt:

echo %GDAL_DATA%

the path to which the variable points to is: C:\Users\me\Miniconda3\Library\share\gdal

This folder contain a file named gcs.csv where the coordinate systems info should be so I do not know what is the problem? I found this discussion:

https://github.com/mapbox/rasterio/issues/1539

but it is too complicated for a beginner like me. Also, there seems to be no solution clearly stated there.

I am using Windows 10, Python 3.7, Miniconda 3, rasterio 1.0.28, and gdal 2.4.2


texts: How can I calculate the distance between multiple points in a network?

texts:
I have multiple points on multiple lines in a network and I need to know the length between these points along the network. I need an OD matrix, like the QNEAT3 distance matrices but for the distance along the network. It's for multiple points, so point and click wouldn't be efficient. Ideally the solution is a similar to the Network Analysis > Shortest Path (point to layer), but it's layer to layer.
Please only suggest QGIS solutions.


texts: Geo-processing Package for ArcGIS Runtime 10.2

texts: I am using ArcGIS 10.4.1 And Java SDK for ArcGIS Runtime 10.2. I created a geo-processing package in ArcMap that is exactly same with the the "SimpleBuffer.gpk" available with the sample dataset. I am using scratchworkspace and set its location according to my output. The sample dataset gpk is running fine and creating tasks in arcgis runtime created server. But when i provide my custom gpk in ArcGIS Runtime it does not generate any tasks.What possible reasons can be for arcgis runtime not creating tasks.  Is there any online tutorial available that would help me in knowning more about creating geo-processing packages.  All of my data is in the same folder to that of gpk. My model builder diagram is exactly same to that of sample buffer dataset.


texts:
I'm trying to import some Landsat 8 GeoTIFFs into GRASS (in QGIS), and I can't work out why I'm getting this error, or how to fix it. No amount of Googling has been able to tell me anything at all. I am using Landsat images as they come, straight from USGS - in other words I haven't modified the images I'm trying to import in any way. The full command output is below.

r.in.gdal input=D:_Eden\landsat\LC80910842015057LGN00_B2.TIF output=b2_GRASS -o

￼Over-riding projection check

￼Fixing subtle input data rounding error of south boundary (-3.94943e+006>1.26566e-007)

￼G_set_window(): North must be north of South

Finished with error

I'm running QGIS 2.10.1-Pisa with Grass 6 v2, and I installed using the Windows Standalone installer. I get the same error when trying to import the same files into the 'standalone' GRASS GIS 6 application.

This is the output I get from g.region -p:

projection: 3 (Latitude-Longitude)
zone:       0
datum:      wgs84
ellipsoid:  wgs84
north:      34:23:23.64S
south:      34:34:23.52S
west:       148:39:18E
east:       148:54:14.4E
nsres:      0:00:00.896576
ewres:      0:00:00.8964
rows:       736
cols:       1000
cells:      736000

Incidentally, it's a southern hemisphere location, as you'd probably guess from the above.

Does anyone have any ideas?


texts:
I'm a new python user trying to build a script tool for an ArcPro project.  The Tool will allow a user to select specific field data from a feature class and then export it as a new feature class to a new .fgdb.  Here are the exact final deliverables:
The tool should contain five parameters (as shown in the image on the right): special notes for parameters setup are shown in paranthesis

input (Data Type: Feature Layer, Type: required)
field (Data Type: Field, Type: required, Dependency: input)
a query expression (Data Type: SQL Expression, Type: optional, Dependency: input)
output geodatabase (Data Type: Workspace, Type: required, Filter: Local Database)
overwrite checkbox (Data Type: Boolean, Type: optional)

The tool should first ask for a input feature layer/feature class.
Then, ask a field from the feature class (set the dependency of the parameter to be the input features)
The tool should also ask for an optional query expression. if specified by user, the tool will export only features that satisfy the query. if not specified by user, use all features in the feature class.
The tool should then export the features as a single feature class for each unique value in the specified field parameter.
The tool should also include a overwrite checkbox.
Each unique value should get a corresponding unique name for the feature class associated with the value. For simplicity, just use the input feature class's name plus a unique id.

My problem seems to stem in my query and final arcpy.SelectLayerByAttribute_management()function.

    Here is my code:
    ##Start by importing the needed drivers
    import arcpy
    import numpy as np

    #Next set up a workspace:
    gdb_worksp = r"SomeLocationOnComputer"
    arcpy.env.workspace = gdb_worksp

    #FeatureClass Parameter
    #Input Feature Class parameter
    input_fc = arcpy.GetParameterAsText(0) #The user's featureclass that they want to use.

    #Field parameter
    input_field = arcpy.GetParameterAsText(1) #The field in the user's featureclass.
    field_arr = arcpy.da.FeatureClassToNumPyArray(input_fc, [input_field])#Converts into NumPyArray
    print("Now the array is not 'structured'.")
    user_field_arr = field_arr[input_field]  #now the array is not 'structured'
    #print(user_field_arr)

    uni_val_arr = np.unique(user_field_arr) #Now it only shows the unique values in the array.
    #print(uni_val_arr)

    #SQL Query parameter
    query_exp = arcpy.GetParameterAsText(2)

    #output_gdb parameter
    output_gdb = arcpy.GetParameterAsText(3)

    #Overwrite parameter
    overwrite = arcpy.GetParameter(4)

    if overwrite is True:
       arcpy.env.overwriteOutput = True
    arcpy.AddMessage(str(overwrite))

    input_lyr = arcpy.SelectLayerByAttribute_management(uni_val_arr, "NEW_SELECTION", query_exp)

      i = 0  #This is the counter
      for value in uni_val_arr:
         query = """"{}" = '{}'""".format(input_field, value)
         input_fc_name = arcpy.Describe(input_lyr).name
         output_name = input_fc_name + "_" + str(i)
         output_fc = output_gdb + "\\" + output_name
         arcpy.Select_analysis(input_lyr, output_fc, query)
         i += 1

'''
If anyone can offer any suggestions I would really appreciate it.


texts: What projection to use to get a game map into QGIS?

texts:
I have built and maintain a SWTOR game map atlas. At the moment, it's all static jpg maps, but I would like to get the info into an actual GIS software to get more functionality and make my website more mobile-friendly.

I assume the first step is to get the game map into a GIS software. I have a png of the map, and I have game map coordinates for a few points on that map. I used the QGIS georeferencer plugin; however, no matter which projection I choose (I tried UTM 19N, UTM 1N, custom eqc I found in this post, even EPSG 2046), my map always shows upside down and/or inverted left-right.

I've read around the internet quite a bit, but still don't understand what I'm doing wrong.

What projection should I be using? Is there a way to define a custom projection and control which axis is pointing upwards (I tried +axis=end without luck)?

For extra clarity, the map:

And the result after attempting to georeference with UTM 1N (same effect with ALL projections I tried):


texts: Google Earth Engine: Layer error: Description length exceeds maximum

texts:
I've imported two shape files into Google Earth engine and am trying to combine them to form a feature collection. I do this using the code;

var ROI = ee.FeatureCollection([
ee.Feature(OSPAR_II.geometry(), {'id': 'OSPAR_II'}),
ee.Feature(OSPAR_III. geometry(), {'id': 'OSPAR_III'})
]);

This seems to print fine, but when I try and add this layer to the map;

Map.addLayer(ROI, {}, 'ROI');

I get the error ROI: Layer error: Description length exceeds maximum.

What is causing this? I've done something similar before and never had this issue.

My code can be found here


texts: batch buffer with different distances with python 2.7 in arcgis 10.1

texts:
I am trying to write a "smart" way in Python to process multiple buffer to lines, polygons and points with different distances. I use Python and ArcMap 10.1.
I can figure out to separate lines, polygons and points in different Feature Datasets and then to apply 3 batches with different "side type" for each Feature Dataset, but I cannot figure out how to apply different distances to each feature class contained in each Feature Dataset.
I saw the post "batch buffer with python in arcgis 10 which is somehow helpful until a certain extent, but I need a further step for the different distances.

Any suggestion out there?

Best regards


texts:
I am having a problem with a contour file I got from a municipality. Two, actually.
I was using 3.16 and it would just hang up at 99% and then fail. This happened with both files.
I tried it again with 3.18 with one file and got this error. But I had chosen structured line, and was using a no-spaces file path on my local desktop, unlike a previous poster. (See images)

I have 16GB of RAM on a Windows 10 platform.
I can't create a DEM for the life of me, and it's stopped my work entirely. Any ideas?


texts:
The new iLife has map functionality, what's the source of the map data used?

Further pictures on engadget, which show a city-level view.


texts:
My pgRouting shortest path algorithm is able to take restrictions and oneway streets into account. Threrefore, OSM data was imported into a PostGIS database using the tool pypgroutingloader. Sadly, the left turn at major intersections is still not accurate. Using a different import-tool (like osm2pgrouting) does not help either.
This is, what a correct route (south to west) looks like (using OSRM).

My setup currently generates the shortest path like this:

As you can see, it is not violating the oneway-rule. It illegaly turns left at the node after the driveway.
How could I get rid of that behaviour? I already thought about calculating the angle between segments and ristricting sharp turns... But I have no idea about how to do that.


texts:
I have been trying for hours to add a color scale to my vector layer in order to avoid showing all those boxes in my legend.

But I don't find any solution. I tried the Geoproc plugin without success, it only works with rasters and my layer is a vector. I also tried generating .qml files but the Geoproc plugin did not accept them.

Does anyone have an idea?

I would like a color scale like this but in red:


texts: ol.style.Rule is implemented in v3.0.0-beta.1 and not in v3.0.0

texts: Can we use ol.style.Rule and ol.expr.parse in  OpenLayers-v3.0.0? I cannot find them but they are implemented in  OpenLayers-v3.0.0-beta.1.


texts:
I'm very new in using GDAL, and trying the tutorial given here.

But with a piece of code with only

GDALRegister_GTiff();
GDALDataset  *poDataset = (GDALDataset *) GDALOpen(argv[1], GA_ReadOnly );

if( poDataset != NULL )
    GDALClose( (GDALDatasetH) poDataset );

Valgrind says there is a memory leak with 32 bytes still reachable. Is there a memory leak in the GDALOpen function?

Edit: Valgrind's output

==25399== 32 bytes in 1 blocks are still reachable in loss record 1 of 1
==25399==    at 0x4C2CC70: calloc (in /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
==25399==    by 0x660B68F: _dlerror_run (dlerror.c:141)
==25399==    by 0x660B197: dlsym (dlsym.c:70)
==25399==    by 0x51A7F27: GTiffOneTimeInit (geotiff.cpp:10764)
==25399==    by 0x51B5382: GTiffDataset::Open(GDALOpenInfo*) (geotiff.cpp:6198)
==25399==    by 0x53C17E7: GDALOpenInternal(GDALOpenInfo&, char const* const*) (gdaldataset.cpp:2318)
==25399==    by 0x53C1A17: GDALOpenInternal(char const*, GDALAccess, char const* const*) (gdaldataset.cpp:2267)
==25399==    by 0x4010D7: main (main.cpp:18)

Line 18 is:

GDALDataset  *poDataset = (GDALDataset *) GDALOpen(argv[1], GA_ReadOnly );

However, I've followed the tutorial, and my program ends with:

/* Once we're done, close properly the dataset */
if( poDataset != NULL )
    GDALClose( (GDALDatasetH) poDataset );

I just observed that if I give a filename directly as a string in GDALOpen instead of using a command line parameter, I don't have the leak anymore.


texts: Split a tmap legend and display in two locations

texts:
I would like to split the legend of a tmap object, so that one of the two elements is in the top left, and one in the top right. I have been through the documentation of tm_legend, but could not find any hint towards that.

The only thing I can imagine is to plot the two legends separately manually, but it would imply to always know the breaks used for each element and it would not be automatic.

Another solution would be to not split the legend and display it on the side, but I would like to gain some space on the figure and display the legend where there are gaps in the maps displayed.

I have copied a MWE from this webpage.

library(tmap)

data(World, metro, land)

tmap_mode("plot")
## tmap mode set to plotting
tm_shape(land) +
  tm_raster("elevation", palette = terrain.colors(10)) +
  tm_shape(World) +
  tm_borders("white", lwd = .5) +
  tm_text("iso_a3", size = "AREA") +
  tm_shape(metro) +
  tm_symbols(col = "red", size = "pop2020", scale = .5)

What would be an efficient way of plotting the two legend elements separately?


texts: How to save WMS getFeatureInfo points to a layer?

texts:
This may be a poorly worded question, because I am not sure of the exact terminology, but I'll give it a shot.

I have imported a vector layer from a WMS server into QGIS. I am able to query the objects on this layer with the Identify tool and get (most of the time) information about that object.

I have another (local) layer of points that I've created. I would like to automate that query process such that a new attribute is created in the point layer from a query to the WMS at that point.

I can't find any information on this, but if the data is made available to the Identify tool my logic tells me that this should be possible. I think I am lacking the appropriate search terms. Any help is appreciated.

Thank you!


texts: Convert lat/lon to ZIP/Postal code using Python


texts: OpenLayers + GeoServer + Oracle: WFS-Transaction update error

texts:
I am trying to implement wfs-t in a map application. tried using ArcGIS server wfs-t, and didn't make it since there is a lock mechanism while trying to edit a feature.ArcGIS WFS-T + OpenLayers; then decided to move on using GeoServer... unfortunately, stuck again. sincerely wish somebody in this site can give me some advice. thanks a lot!

what I did in JS code is:

update the value of an attribute, and set the state of that feature as OpenLayers.State.Update;

feat.attributes['MATERIAL'] = 'CMPtest';
feat.data['MATERIAL'] = 'CMPtest';
feat.state = OpenLayers.State.UPDATE;

then:

using SaveStrategy.save();

the error response text is:

 <?xml version="1.0" ?>
<wfs:WFS_TransactionResponse version="1.0.0" xmlns:ogc="http://www.opengis.net/ogc" xmlns:wfs="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://XXXX/geoserver/schemas/wfs/1.0.0/WFS-transaction.xsd">
  <wfs:InsertResult>
    <ogc:FeatureId fid="none"/>
  </wfs:InsertResult>

  <wfs:TransactionResult>

    <wfs:Status>

      <wfs:FAILED/>

    </wfs:Status>

    <wfs:Message>
      update error
    </wfs:Message>

  </wfs:TransactionResult>
</wfs:WFS_TransactionResponse>

the request xml is:(from Chrome console)

<?xml version="1.0" ?>
    <wfs:Transaction service="WFS" version="1.0.0" xmlns:wfs="http://www.opengis.net/wfs" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.0.0/WFS-transaction.xsd">
      <wfs:Update typeName="feature:SAMPLE_FEATURE" xmlns:feature="http://www.opengeospatial.net/cite">
        <wfs:Property>
          <wfs:Name>
            SHAPE
          </wfs:Name>
          <wfs:Value>
            <gml:Point srsName="EPSG:26912" xmlns:gml="http://www.opengis.net/gml">
              <gml:coordinates cs="," decimal="." ts="">
                401550.55,4647703.82
              </gml:coordinates>
            </gml:Point>
          </wfs:Value>
        </wfs:Property>
        <wfs:Property>
          <wfs:Name>
            MATERIAL
          </wfs:Name>
          <wfs:Value>
            CMPtest
          </wfs:Value>
        </wfs:Property>
        <wfs:Property>
          <wfs:Name>
            CENTER_LON
          </wfs:Name>
          <wfs:Value>
            -112.188498069
          </wfs:Value>
        </wfs:Property>
        <wfs:Property>
          <wfs:Name>
            CENTER_LAT
          </wfs:Name>
          <wfs:Value>
            41.9749785934
          </wfs:Value>
        </wfs:Property>
        <wfs:Property>
          <wfs:Name>
            CENTER_Z
          </wfs:Name>
          <wfs:Value>
            1353.59340947
          </wfs:Value>
        </wfs:Property>
        <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
          <ogc:FeatureId fid="SAMPLE_FEATURE.1"/>
        </ogc:Filter>
      </wfs:Update>
    </wfs:Transaction>


texts: SelectLayerByLocation keeps giving error 000840: The value is not Boolean

texts:
I've seen a lot of questions asked and answered about trouble with the Select Layer By Location tool. Most ask about invalid input errors and the answer generally is to make sure to use layers (as apposed to feature classes). My problem is that I keep getting an error that some value is not Boolean, which I find strange since none of the inputs require boolean values.

I have two input files. Input1 is a shapefile of building footprints. Input2 is a shapefile of various shaped polygons. What I want to do is iterate through the polygons in Input2 and use each polygon to select the building footprint polygons in Input1 to do some work on them. I can't get past the Select By Location tool error to complete the work. Here is my code example:

with arcpy.da.SearchCursor("Input2","ID") as cursor:
  for row in cursor:
    value = str(row[0])
    field = "ID"
    exp = field + "=" + value
    arcpy.SelectLayerByAttribute_management("Input2","NEW_SELECTION",exp)
    arcpy.SelectLayerByLocation_management("Input1","WITHIN","Input2")
    ...do work on selected building polygons

I have the two inputs active in my arcpy map so they should be interpreted as layers. This is proven true since the Select Layer By Attribute tool works. I put the code into the python window, run it, and get this error message:

...in SelectLayerByLocation Execute Error: Failed to execute. Parameters are not valid. ERROR 000840: The value is not Boolean. Failed to execute (SelectLayerByLocation).

I gave the beginning of the script for context, but I've tried running the Select Layer By Location tool by itself in the python window and still get that message. When I use the Select Layer By Location tool from the tool box, I have no problem and it works. I've tried right clicking on the result, Copying As Python Snippet, pasting into the python window, and running, but I get the boolean error. I've tried using the Make Feature Layer tool in the python script (even though I shouldn't have to) with the same result. I've tried restarting.


texts: Pivot Table deleting field values when pivoting table on single field in ArcMap?

texts:
I have a polygon shapefile (multipart) that I would like to pivot on a single field using ArcMap 10.4.

Everything works well but the final table miss the values of the pivoted field (that are important as well as others).

Can anyone explain me why it happens and how to solve it?

I checked "GRIDCODE" and "CODE" with "CODE" as pivot

Checking "CODE" as Input and "GRIDCODE" as Pivot:


texts: ArcGIS Near tool returning incorrect results?

texts:
I'm having some odd issues with the Near tool in ArcMap 10.4.1 using Geodesic method, no location or angle specified. The tool works, but returns near FIDs for polygons that are absolutely not the nearest. For instance, I can visually see that polygons from the input feature and near feature are adjacent, but the NEAR_FID is for a different polygon.

I've repaired geometry on all involved parties, converted multi-part to single-part for the near features, and tried a variety of search radius sizes.


texts:
I'm trying to clip a raster using a feature set. I'd like to then have a new value added to a column in attribute table of the feature set based on the majority of the pixels in the raster.

In this illustration, raster values could be corn (y), soy (g), no data (trans). I'd like to clip with the feature and then have the "crop" column in the clip feature attribute table updated to be "soy" if it's the majority value.

QGIS solutions? ArcGIS?


texts: Overpass QL for Polygons of Land between streets

texts:
I am trying to create a Overpass Query that would give me polygons of Land that is NOT a highway/street/waterbody/(anything place that cannot be residential).
The area Im looking does not have much mapping data. I wanted to have all residential areas and buildings but they arent mapped yet.

I need an alternative that would be a good guestimate of where these residential areas might be. (Not roads/waterbodies/parks)

I just can't seem to polygons of the spaces between the streets in the picture below.

UPDATE

Link to Overpass Link http://overpass-turbo.eu/s/j3I
This will fetch the roads and water bodies, I don't want the land, I just want the spaces between those roads


texts:
I am creating a plugin but my nested loop does not work well and I do not know why. How can I do this?
In PyQGIS the loops are different from the Python?
I have 2 inner loops inside a for loop, but my console ignores something, and I don't understand why?
In this case the loop ignores the points that I have on TFEA iterator.
I wrote the following code:
for f in intersecting_features:
    print("a")
    layer.select(f.id())
    TFEA = point_of_study_feature.getFeatures()
    dist = 500
    for p in TFEA:
        print("b")
        analyze_feature_point = processing.run("qgis:pointsalonglines",{'INPUT':processing.run("native:saveselectedfeatures", {'INPUT':layer,'OUTPUT':'memory:'})["OUTPUT"],'DISTANCE':1,'START_OFFSET':0,'END_OFFSET':0,'OUTPUT':'memory:'})["OUTPUT"].getFeatures()

        for q in analyze_feature_point:
            dist_data = distance.measureLine(p.geometry().asPoint(), q.geometry().asPoint())
            if(dist_data < dist):

                dist = dist_data
    if (f["Z_MAX_ED"] - ((dist)*0.1)) >= select_feature["Z_MAX_ED"]:

        f.setAttributes([f["IDP"], dist, f'Sombreador: {f["IDP"]} (Z_MAX_OBST= {f["Z_MAX_ED"]}, dista {dist}m); Z_MAX_SOMBREAMENTO={(f["Z_MAX_ED"] - ((dist)*0.1))}m >= {select_feature["Z_MAX_ED"]}m.'])
        pr.addFeature(f)


texts: How do I get a Python list of field names in PyQGIS 3.x?


texts: Inserting values in a table depending on spatial relation between objects in to tables

texts:
I have a new and empty column in a table and I want to insert values, using postgis depending on the spatial relation between the geometries in two tables.
I have one table (points) with benches and another table with cadastres. I now want to add a 1 (one) in the benches table where objects in the benches table intersects cadastres that public owned and a 0 (zero) where they intersect cadastres owned by anybody else.
The cadastres table is EPSG: 25832
The benches layer is EPSG: 4326
I have tried the following:
SET public =CASE WHEN ST_Intersects(a.the_geom,b.the_geom) THEN 1
ELSE 0 END
FROM schema.object_bench a
INNER JOIN (SELECT * FROM schema.cadastres  WHERE owner = 'public') b
ON ST_intersects(ST_TRANSFORM(a.the_geom,25832),b.the_geom);

When I run the query there is no syntax error message but no rows are updated, which not should be the case since all the objects in the benches in table intersects objects in the cadastres table.
Is it the ST_Transform operation that fails or is the query in general?


texts: What projection to use to create aligned GB and Ireland OS Grid map

texts:
I would like to create a map like the following:

It displays the 100km grid square systems used by both UK and Ireland.

I have created the 100km grid squares using the TomBio tools plug in, however, the issue I'm having is with finding a projection over both the UK and Ireland that ensures the grid squares align up as in the image above with neither Ireland of UK rotated such as in the image below, which happens when using either 27700 or 29902. Is this possible?


texts:
I am using GeoServer and seed tiles on my server. The tiles are created successfully but I don't know which pattern the directory structure is following... (i.e. .../EPSG_4326_05/0_0/00_06.png).

I want to use the tiles in a OpenLayers application and there I want to use a OSM source which is using the XYZ-pattern which is commonly used as URL pattern for tile-serving. Is there a way to tell the GeoServer it should create the tiles with the XYZ structure?


texts: Identifying attributes from related tables

texts:
I've created a relationship class and created a 1-m relationship between my feature class and table. When identifying the results in ArcGIS 10.2 it returns the details using the "OBJECTID" field. Is it possible to query the results and show the "FlightID" rather than "ObjectID"?

In the example provided in the following, the 7266 ("OBJECTID") would be replaced by 16-284 ("FlightID"):


texts: Why does exporting labels with buffers take so long?

texts:
I set white borders for text (city names) in qgis. If I export the map without white borders it exports quite quickly (png 600 dpi). If I set the white borders it takes hours! Is there any way I can avoid this?

Best

the qgis Version is Lisboa. Operation system Linux. The text labels are set in the regular Text Tab of the Layer-Properties. See screenshot http://www.projektidee.org/images/extern/qgis2.png. In the map composer I didn't change any settings.


texts: Is there any way to use callout windows or boxes in MapInfo with lead lines to label a point layer?


texts:
QGIS 3.8.3 running in Ubuntu 18.04.3 LTS

I've installed Whitebox Tools and the Whitebox plugin and am able to obtain, for example, data on a .las file using 'Lidarinfo', so the system appears to be working (at least at this level).

However, attempting to use LidarNearestNeighbourGridding produces the following:

  Algorithm 'LidarNearestNeighbourGridding' starting… Input parameters:
  { 'exclude_cls' : '1,3,4', 'input' :
  '/home/nick/Desktop/SU5771_P_9194_20131130_20131201.las', 'maxz' :
  None, 'minz' : None, 'output' :
  '/home/nick/Desktop/SU5771_P_9194_20131130_20131201.tif', 'parameter'
  : 0, 'radius' : 2.5, 'resolution' : 1, 'returns' : 0 }

  WhiteBoxTools command: "/home/nick/whitebox_tools/WBT/whitebox_tools"
  --run=LidarNearestNeighbourGridding --input="/home/nick/Desktop/SU5771_P_9194_20131130_20131201.las" --parameter="elevation" --returns="all" --resolution=1.0 --radius=2.5 --exclude_cls="1,3,4" --minz=0.0 --maxz=0.0 --output="/home/nick/Desktop/SU5771_P_9194_20131130_20131201.tif" -v WhiteBoxTools command output:

  Welcome to LidarNearestNeighbourGridding *

  Performing interpolation...

  Reading input LAS file...

  No points found in
  /home/nick/Desktop/SU5771_P_9194_20131130_20131201.las

  Finished interpolating homenickDesktopSU5771_P_9194_20131130_20131201
  (1 of 1)

  Elapsed Time (including I/O): 0.146s

  Execution completed in 0.29 seconds Results: {'output':
  }

  Loading resulting layers The following layers were not correctly
  generated./home/nick/Desktop/SU5771_P_9194_20131130_20131201.tifYou
  can check the 'Log Messages Panel' in QGIS main window to find more
  information about the execution of the algorithm.

"No points found in /home/nick/Desktop/SU5771_P_9194_20131130_20131201.las" would appear to be the operative part of this.

The file is a UK Environment Agency .laz file unzipped using laszip and was supplied already classified.

I'm very new to processing lidar data but I have produced DTMs using LAStools and also using PDAL.

Am I missing something here?


texts: How to connect to ArcGIS Server on the other machine?

texts:
There is the ArcGis Server on the machine 192.168.**.22

Machine name is MAPSERVER.

User name and password for ArcgIS Server are "test" and "testtest" (for example)

ArcGis Server options are: clusters - name "default", computer "MAPSEVER", protocol "TCP port 4004".

I want to connect to this ArcGis Server. I used code

ESRI.ArcGIS.ADF.Identity identity = new ESRI.ArcGIS.ADF.Identity("test", "testtest", "default");
ESRI.ArcGIS.ADF.Connection.AGS.AGSServerConnection agsconnection;
agsconnection = new ESRI.ArcGIS.ADF.Connection.AGS.AGSServerConnection("MAPSERVER", identity);
agsconnection.Connect();
IServerObjectManager SOM = agsconnection.ServerObjectManager;

I get error `

  "Identity error"


texts:
I'm building a plugin from QGIS, and I need to extract the values from the quartiles from a band in a raster layer.
Does anyone knows how it is done?


texts:
What I'm trying to do is find a road (or a portions of road) from one country that is similar to another road(or portion of a road) from a different country and the road property I'm comparing is relative changes in elevation.

For example, if I have a 100km road stretch in Kenya, looking at how the elevation changes every 1km, I would like to know if there is a either a continuous 100k road stretch in Namibia with very similar elevation changes or a smaller stretch (if the full length doesn't exist).

Since I have no prior GIS experience, I started using the Google Maps API in Javascript (and it's Elevation API), but I'm open to any suggestion.

At the moment I'm considering taking a possibly very naive approach:

store input road data from country #1 as a 1 dimensional list of delta elevations (every 1km)
index roads from country #2 as 1D lists of delta elevations
(somehow) search for input road in the indexed road data. Either iterated through the indexed 2 lists and sort by the average difference between two sections or attempt 1D template matching ?

Is what I'm attempting achievable ?

Am I asking the right question (or I should look at comparing roads in a different way) ?

If so, is there an efficient/elegant solution to searching similar road areas based on relative elevation changes ?


texts: returning extent for the selected features from microsoft sql server in geoserver

texts:
I have data in the GeoServer coming from Microsoft's SQL server. I want to use CQL filter with WMS to query and get 4000k records as a result. I need to know is there any way my CQl filter with WMS or WFS or WPS can return extent for the selected 4k records in the resultant response?If the result is 2k records, then it should return the extent for those.For ex,

http://localhost:8080/geoserver/wfs?service=wfs&version=1.1.0&request=getfeature&typename=topp:states&CQL_FILTER=STATE_NAME=%27Idaho%27&resulttype=hits

returns 1 record as below in XML

<wfs:FeatureCollection xmlns:it.geosolutions="http://www.geo-solutions.it"    xmlns:cite="http://www.opengeospatial.net/cite" xmlns:ogc="http://www.opengis.net/ogc" xmlns:tiger="http://www.census.gov" xmlns:sde="http://geoserver.sf.net" xmlns:topp="http://www.openplans.org/topp" xmlns:wfs="http://www.opengis.net/wfs" xmlns:ows="http://www.opengis.net/ows" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:KFRI="http://www.kfri.res.in/" xmlns:sf="http://www.openplans.org/spearfish" xmlns:gml="http://www.opengis.net/gml" xmlns:nurc="http://www.nurc.nato.int" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" numberOfFeatures="1" timeStamp="2016-08-01T11:45:42.443Z" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"/>

In addition,i need the extent also for the feature in the response. And please note that if the response returns 4 features also then it should return one extent for all those 4 features instead of returning extent for each feature.


texts:
I am using Leaflet with OpenLayers and I have the following code:
The all code I used is the following :
OSMaps = tileLayer("https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png", {
  attribution:
    'Map data &copy; OpenStreetMap contributors, CC-BY-SA, Imagery © CloudMade',
  maxZoom: 50,
  detectRetina: true,
});

wMaps = tileLayer(
  "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}",
  {
    maxZoom: 30,
    detectRetina: true,
    attribution:
      "Tiles &copy; Esri &mdash; Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community",
  }
);
dMaps = tileLayer(
  "https://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png",
  {
    attribution:
      '&copy; OpenStreetMap contributors &copy; CARTO',
    subdomains: "abcd",
    maxZoom: 30,
  }
);

// Marker center of France
center = marker([46.4547, 2.2529], {
  icon: icon({
    iconSize: [25, 41],
    iconAnchor: [13, 41],
    iconUrl: require("../../../../../../node_modules/leaflet/dist/images/marker-icon.png"),
    shadowUrl: require("../../../../../../node_modules/leaflet/dist/images/marker-shadow.png"),
  }),
});

// Layers
devicesLayer = new LayerGroup();
options = {
  layers: [this.OSMaps],
  //zoom: 50
  minZoom: 0,
  maxZoom: 50,
};

layersControl = {
  baseLayers: {
    "Wikimedia Maps": this.wMaps,
    "Night Maps": this.dMaps,
    "Street Maps": this.OSMaps,
  },
  overlays: {},
};

The way of adding points is the following :
device_marker = marker(deviceLatLng, {
  icon: icon({
    iconSize: [32, 46],
    iconAnchor: [22, 46],
    iconUrl,
  }),
  //title: item.name
  //Aqui de pone el titulo
});

this.markersMap.push(device_marker);
this.map.addLayer(this.markerMap);

When I try to set Leaflet in a map and I set zoom to the map I obtain the following error:
Failed to load resource: the server responded with a status of 404 ()


texts:
When I'm loading a CSV into QGIS, I get this error code, see below.
Does anyone know what this means?
This also is an issue when I attempt to run NNJoin Plugin too.
An error has occurred while executing Python code:
  NameError: name 'QMessageBox' is not defined
  Traceback (most recent call last):
    File "/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/NNJoin/NNJoin_gui.py",
      line 367, in layerlistchanged
      QMessageBox.information(None, NameError: name 'QMessageBox' is not defined

Python version: 3.7.7 (default, Sep 22 2020, 10:25:18) [Clang 12.0.0 (clang-1200.0.32.2)]
QGIS version: 3.16.1-Hannover Hannover, 3ecbe06f78

Python Path:
/Applications/QGIS.app/Contents/MacOS/../Resources/python
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins
/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/Rtree-0.9.4-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/cftime-1.2.1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/snuggs-1.4.7-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/Pillow-7.2.0-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/lib-dynload
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/rasterio-1.1.5-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/matplotlib-3.3.0-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/netCDF4-1.5.3-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python37.zip
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/affine-2.3.0-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/patsy-0.5.1-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/scipy-1.5.1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/numba-0.50.1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/geopandas-0.8.1-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/pandas-1.1.0-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/Fiona-1.8.13.post1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/statsmodels-0.11.1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/opencv_contrib_python-4.3.0.36-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/click_plugins-1.1.1-py3.7.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/numpy-1.19.1-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/GDAL-3.1.2-py3.7-macosx-10.13.0-x86_64.egg
/Applications/QGIS.app/Contents/MacOS/lib/python3.7/site-packages/pyproj-2.6.0-py3.7-macosx-10.13.0-x86_64.egg
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/DigitizingTools/tools
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/DigitizingTools
/Users/joanlee/Library/Application Support/QGIS/QGIS3/profiles/default/python/plugins/mmqgis/forms
/Users/joanlee/Camphora/15. ITTC/Maps


texts: Where can I get seamless 500dpi DRGs for Oahu (Honolulu County)?

texts:
USGS makes 500dpi DRGs but doesn't seem to offer them for download anywhere. There are several vendors that I can find offering seamless DRGs but only in 250dpi. Where can I find 500dpi seamless DRGs for Oahu?

GeoPDFs will not work-- I have already tried importing them to Photoshop, exporting high res TIFFs, and georeferencing those, but I really need the edges to match perfectly.


texts: Digitizing contour lines from paper maps using QGIS?

texts: Is it possible to give digitized altitude lines (from paper map) heights using QGIS?


texts: Need to assemble a matrix from several ASCII GRD files

texts:
I am totally new to GIS software.

I would like to consolidate 205 ASCII GRDs files (the file format of Golden Software's "Surfer") into a single big GRD file or just a big matrix of (about 2400 rows by 2400 columns) I can manipulate within Excel.

Each of the 205 GRDs contain mean windspeed (32 bit float) information of 150x150 cells. I could import these GRDs one by one in Excel but I thought there could be a faster way to do this (I actually have 5 other parameters in this same GRD format).

I was able to put these 205 GRD files in one Arcmap layer using the "Mosaic to New Layer" tool but I didn't find a way to export this layer to a single GRD. I haven't tried the same in QGIS.

Any ideas ?


texts: Can anyone suggest best Suitable Algorithm for Digital Elevation Model Generation? Currently I am working on Literature survey of DEM creation so i need to show comparison of 3-4 Algorithmic Methods used to create DEM & Best suitable method among them.


texts: Batch processing Cell Statistics in ArcGIS using Python - ERROR 000864

texts:
The Below script is designed to take all raster (.txt ascii rasters) in a specified folder, add them to a list with the file paths, use the file paths stored in the list in the cell statistics tool to get the mean. The script then iterates on to the next folder until the while function is satisfied.
There is a folder for every year from 1969 - 2014.

This script has previously ran successfully, and i am sure i haven't changed anything since i last ran it but today when i ran it i got the error:

"ERROR 000864: Input raster's or constant values: The input is not within the defined domain."

This error usually is regarding to complications with floating bit and integer but i thought that Cell Statistics could handle both types - i have also completed cell statistics on these raster's manually so i know they can be used.

Arcpy indicates the error is on line 23 which is the cell statistics part of the code.

import arcpy
import numpy
import os

file = 1969 #Set the starting year as a variable

while file <=2014:  # set the loop to run up to the year 2014

    endPath = str(file)#get the starting eyar as a string

    #get the contents of the folderstructure in a list
    dirListing = os.listdir('C:\\Users\\Desktop\\MetOffice\\MetOffice Data\\'+ endPath)

    file_Path = 'C:\\Users\\Desktop\\MetOffice\\MetOffice Data\\'+ endPath + '\\'
    file_Path_ext = [file_Path + x for x in dirListing]
    #get the full file path for each of the files and join it to the files.

    #Output Location for the processed shapefiles
    output = "C:\\Users\\Desktop\\MetOffice\\ProcessedData\\"+endPath+".tif"

   # Process: Cell Statistics
   arcpy.gp.CellStatistics_sa(dirListing, output, "MEAN", "DATA")
   file = file + 1 #add one to the year to start the next loop and move
   through the folders
   print(endPath + " Completed" + ", " + str(2015 - file) +" Remaining")#print a update

print("All Years Completed")#tell the user when the task is completed


texts: How do I stop FME from using the first row as the column name when reading XLS?

texts:
I have a number of Excel XLS-files that I want to process with FME. They're each structured in two columns, column A with area names and column B with values for each area.

When reading the file, FME automatically uses the top row as column names. However, this row is actually data as well rather than column headers.

Now I can't figure out how to stop this behaviour in FME. So, how do I expose the columns as generic (A and B for example), and prevent FME from using the top row as column names?


texts: Difference between LAStools, libLAS and PDAL?

texts:
After many readings, I'm still confused about the difference of use and the aims of the following tools for LiDAR data processing and manipulation: LAStools, libLAS, and PDAL.

My understanding so far is that PDAL and libLAS have the same overall goals, and that PDAL is slowly replacing libLAS with enhanced performance and more capabilities. I understand that PDAL is a library that is (and want to be) completely independent from LAStools.

The relation between LAStools and libLAS is the most confusing to me. Are they two completely different separate tools? Does LAStools run on top of libLAS?

Does any of these 3 tools/libraries depends on one another?

Background: I have experience with using LAStools, and have played around with PDAL in the past few weeks.

Edit: From a link comment, I learnt that there is also something called LASlib which seems to be something LAStools runs on. I wish this question can gather a comprehensive overview of these tools for people starting off using LiDAR data and being confused what to go with.


texts:
In my workflow I typically, in Microstation V8 2004 Edition, manually place a shape element around a series of imported coordinates to use as a fence for loading point cloud data in TerraScan.  I am trying to incorporate QGIS into my workflow doing the following:

Importing and buffering my coordinate file in QGIS
Editing the file until my buffer is one continuous polygon
Simplifying the geometries so I am under approx. 100 nodes
Converting the polygon to a line shapefile
Exporting this shapefile to .dxf format for use in Microstation

My issue is that when I attach the .dxf file to my working DGN, and merge it into the master (effectively making it an element), Microstation is not recognizing it as a valid fence and won't constrain the point import to with my desired area.

I am wondering if there is an issue with the geometry, from the way I created the line file in QGIS?

DXF files containing line polygons converted in Global Mapper (from different source data) work as expected - that is, they constrain the data in the same way an that an element does when I create it manually within Microstation. So I have also tried importing my QGIS line shapefile into Global Mapper and exporting it as a .dxf, but this file does not work either - leading me back to it being a QGIS issue.

Anyone have any thoughts or tried to do this before?

I am a relatively new QGIS user, with Version 2.4.0-Chugiak


texts: Does QGIS Server come with time dimension support (WMS, WCS)?

texts:
I wanted to know if QGIS Server allows publishing data with time dimension support, in the case of services such as WMS or WCS?. If this is possible, how can I enable this feature?.

Thanks for your answers.

regards.


texts: QGIS 2.18 Scale not being set up correctly

texts:
I'm working on adding scale based visibility on different layer. I'm currently using :

canvas.setDestinationCrs(crs)
crs.createFromUserInput(name)
crs = QgsCoordinateReferenceSystem()

to set the CRS of my project. However, when launching it with QGIS 2.18, it won't load until I go to the project's settings and click OK in the CRS panel. This problem does not apply to QGIS 3
Below is a GIF showing how to solve it using the GUI:


texts: Comparing results of cross K-L functions for point patterns

texts:
Given marked point patterns with types a, b, and c, spatstat in R allows to easily calculate cross K(a,b) and L(a,b) functions. Envelopes can be used to estimate the significance at a given alpha (0.05) for each pair at each distance, for example at 500 meters (a,b) in envelope, (a,c) in envelope, (b,c) not in envelope.

How can I see whether the K/L functions for two pairs ((a,b), (a,c)) are significantly different to each other? This is to state that for example trees of type a and b cluster in a significantly different way to a and c at a given alpha level.


texts:
I'm new to GIS in general and field calculator in particular.

I have a basic attribute table with three columns. I'd like to assign values to column C based on what is in column A. If Column A = 8 then I need column C to read "pine forest". Column C is set for text with length=100 so that should not be a problem.

I had thought I could do this by selecting all of the records where column A = 8 and then using field calculator in column C to assign them the label "pine forest". However this results in an error message saying that the calculated value is invalid because it is either too large or because it is a string, which are both wrong.


texts:
I have photogrammetry DEM with 3 meter resolution. I need grid point with 30 meter resolution with the 3 meter accuracy.

Case-1,
I created the the 30x30 meter grid points and extract the elevation from DEM.

Case-2,
I interpolate the 3x3 meter DEM to 30x30 meter raster by using near neighbor interpolation.

Case-3,
Create TIN from Raster and create raster 30x30 meter.

Could you guide which is the best method to interpolate high resolution to low resolution raster?


texts: PostGIS: Using ST_CoveredBy (or similar) to speed up Difference on large datasets?

texts:
I need to get the Difference of a single large-area polygon (input layer) and a very numerous set of polygons (overlay layer). The overlay layer has ~500k polygons, of which ~200k intersect the input layer. I am trying to use PostGIS to do this in the hope that it can be much faster than e.g. running a Difference algorithm in QGIS (which is unfeasibly slow).
I've got the basic approach to do this with ST_Difference from here, i.e. just run ST_Difference on a left join on ST_Intersects (using Coalesce to avoid losing polygons from the input layer that don't intersect with the overlay layer - although as it happens in this particular case there are no such polygons to lose as there is only a single input polygon which does intersect).
So my current query is:
SELECT COALESCE(ST_Difference(input.geom,overlay.geom),input.geom)
FROM input LEFT JOIN overlay ON ST_Intersects(input.geom,overlay.geom)

Unsurprisingly this query is taking a very long time, as it's calculating ~200k differences.
I've read about speeding up intersections enormously by skipping ST_Intersection where polygons are entirely contained within one another e.g. here and here. I.e. using ST_CoveredBy / ST_Within / ST_Contains etc. (I haven't completely sussed out the difference between these yet.)
I feel like I should be able to do something similar here, as the vast majority of the overlay polygons are entirely contained within the input polygon. But I can't get my head around the query code to use to do the same with ST_Difference.


texts: Get shapefile polygon attribute value at a specific point using Python (e.g. via PyShp)

texts:
I have a shapefile comprising multiple polygons. These polygons have attributes. I am trying to write a script in Python that takes a coordinate point, determines in which polygon it falls, and returns the attributes for that polygon.

So far, using PyShp I've got this, which prints off all the records.

import shapefile
sf = shapefile.Reader("England.shp")
records = sf.records()
print records

So I want to pass in a coordinates for Buckingham Palace, which is inside a polygon which has the attribute 'Area' value of 'London' and have it return 'London'

Can anyone point me in the right direction?


texts: Create a temporary intermediate layer in a processing script

texts:
My algorithm needs a temporary layer in a QGIS processing script. I can use regular pyqgis api to create a temporary layer but want to know if there is a preferred and processing-specific way of doing it. Since processing already creates a temporary directory for the output, can it be used for creating another layer? Can I use processing.core.Vectorwriter  or some other helper method?

Edit:
I use the following right-now and looking for something more processing-specific and preferably writing to processing temporary dir instead of a memory layer.

templayer = QgsVectorLayer('Polygon?crs=%s' % input_layer.crs().authid(), 'temp', 'memory')
dataProvider = templayer.dataProvider()
templayer.startEditing()
dataProvider.addFeatures([features])
templayer.commitChanges()
templayer.updateExtents()