texts: finding a location on a line using a point and a relative distance with arcgis js api

texts:
The problem I'm trying to solve is to find a point on a line given another point on the same line and a distance from the same point. I've been looking on the various RouteTask examples but they mostly tackle the problem of finding a route given two points. Can anyone advise if RouteTask is the right approach for this kind of problem, and how do I do it?

Thanks!


texts: Calculate Percentage population from Vector Grids

texts:
I have a vector grids layer of population (shapefile). Each grid represents population in that grid and have the following attributes:
Area
Population
I want to calculate and display the percentages of population and their respective area of coverage.
For example: which grids (highlight those grids) constitutes the x% of the population (to indicate the area where we should deploy our resources for maximum coverage), how many are they and what is the total area (effectively, number and sum of the attributes of selected grids).
I am trying to explore something that can help in QGIS but haven't found anything yet. Can anyone please recommend something (any tool, script or software as I am not bound in terms of software).


texts: Removing only donut holes from polygons using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.0.

I have a trick for removing the donut holes, which is converting the polygon to line and then creating polygon out of the lines and in the end merging all of them. The drawback of this strategy is that in removes other holes and empty areas which in reality are empty and should be kept empty. See below image:


texts: Projection parameters for Gauss Kruger 7 zone Serbia

texts: I can not set  parameters for Gauss Kruger 7 zone Serbia


texts: QGIS 2.14.3 Processing Dissolve error "local variable 'nElement' referenced before assignment See log for more details"

texts:
When I try to use Dissolve in Processing I get this error:

"local variable 'nElement' referenced before assignment See log for more details"

I have seen this problem reported in:  https://hub.qgis.org/issues/14612

So then I did what is recommended by  Alexander Bruy removed the processing folder from C:\Users\user.qgis2 however I still have the same problem, I did a Full Uninstall (delet all folders...) of Qgis and reinstalled from standalone installer QGIS-OSGeo4W-2.14.3-1-Setup-x86_64 but the problem persist.

I suspect that it is a bug...

PS: when I use the "Dissolve all" check, the process work, its only when I try to dissolve from a selected field that give me the error

QGIS 2.14.3
Windows 10 64bits (1511 10586.318)
Processing 2.12.99


texts: Adding JAXA weather map tiles in Leaflet

texts: Is there any way I can add JAXA rainfall layer as a Leaflet layer. If not what other options are available for rainfall/ppt layer other than OpenWeather map?


texts:
Based on the following code block:

var SphericalMercator = require('sphericalmercator');
var mercator = new SphericalMercator({
    size: 256
})

var bbox = mercator.bbox(
   +req.params.x,
   +req.params.y,
   +req.params.z,
   false,
   '4326');
var test_query = 'SELECT st_asgeojson(wkb_geometry) AS feature FROM ogrgeojson where
ST_INTERSECTS(wkb_geometry, ' + bbox + ')';

I run this as a query and get the following error:

error: function st_intersects(geometry, numeric, numeric, numeric, numeric) does not exist]

hint: 'No function matches the given name and argument types. You might need to add explicit type casts.'

I am not sure how I should be wrapping the bbox. Additionally things that I have tried have resulted in column errors.


texts: Publish multiple tables from different PostgreSQL's schemas on GeoServer

texts:
I have many tables (about 1000) in a postgis which  are seperated in four different schemas. I know that every PostGIS store must be connected with just one schema. Thus, what is the most efficient way to publish these tables from different schemas on GeoServer ?

Is it a good practice to create four PostGIS stores with the same user/password but with different schema or have I to use PostGIS JNDI connection ?
Is there another way to go ?


texts: Running a geoprocessing service on selected data in ArcGIS web appbuilder

texts:
I am looking to run a geoprocessing service on a selected subset of feature data in a web app. I can use the entire feature layer as a parameter, but cannot seem to use only features selected either on the map or in the data table.

The end goal is to select one or more parcels in the web app and run a tool on only the selected features. I have tried using both the feature service and the map service as an input.

Is there a way to do this without resorting to using the interactive drawing input and running a spatial join?


texts: Get all rows that have geometries which intersect with a geometry in any row of a different table in an effective way

texts:
I have two PostGIS tables containing thousands of geographic polygons each (each row has one polygon and numerous other attributes). Both tables represent different data sets but cover the same area of interest.

I'm trying to receive the rows of table A that intersect with rows from table B (meaning that the geometries intersect).

What I tried is the following:

SELECT tableA.polygon FROM tableA, tableB
WHERE ST_Intersects(tableA.polygon, tableB.polygon)=true;

I have no idea if this works or not since the query has not finished since I started it the first time 20 minutes ago. I have to admit that I'm running the query on a machine which is not very powerful.

Is this the correct way to do it? Is there a smarter or more effective way allowing me to reduce the computing time?


texts: Finding furthest Distance from Road using ArcGIS Spatial Analyst?

texts:
I have two vector layers, an administrative boundary layer and a roads layer. They look like the following:

I would like to find the point (within the administrative boundary) that is furthest from a road. I have a spatial analyst licence and I have tried running the Euclidean distance tool but it merely buffers around the coastline.

Do I need to combine both layers first?

Do I need to rasterise them individually and then run the Euclidean distance tool?

I have looked at dozens of other questions on here and I can't seem to find an answer to my problem.


texts: Arcpy: copied python snippet does not work?

texts:
I am trying to convert XYZ-data to a point feature class. I've done this successfully in ArcMap, however, when I copy the python snippet, I get a RuntimeError.

I did make some minor changes to make the code more readable, and to loop through a series of files. Anyway, the same problems arose when I used the original python snippet.

This is my code - for the sake of brevity I have replaced some of the long path names by variables:

import os                    # enables interaction with operating system
import arcpy                 # enables use of tools within the arcpy package
arcpy.CheckOutExtension('3D')

RDnew = "PROJCS['RD_New',GEOGCS['GCS_Amersfoort',DATUM['D_Amersfoort',SPHEROID['Bessel_1841',6377397.155,299.1528128]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Double_Stereographic'],PARAMETER['False_Easting',155000.0],PARAMETER['False_Northing',463000.0],PARAMETER['Central_Meridian',5.38763888888889],PARAMETER['Scale_Factor',0.9999079],PARAMETER['Latitude_Of_Origin',52.15616055555555],UNIT['Meter',1.0]]","#","#","DECIMAL_POINT"
#this is just the coordinate system copied from the python snippet to make it more readible.

for fn in os.listdir('.'):
    if (fn[-4:] == '.pts' or fn[-4:] == '.txt'):
        print fn + ' is processing.'
        print 'Converting XYZ-data to point feature class.'
        pts_shp = fn[0:7] + '.shp'      # using fn[0:7] to get rid of the .pts or .txt extension of the input file
        arcpy.ASCII3DToFeatureClass_3d(fn,'XYZ','shp_points/' + pts_shp,'POINT','1',RDnew,"#","#","DECIMAL_POINT")

This last line is where I get the error.

I had some errors before (among which the infamous Error 999999, which I managed to get rid off. and 'failed to execute function', which has also been solved), but I don't see what else is wrong with my code here.

Could anyone help me out and explain to me what is going wrong?

I already checked for spaces and dots in the paths and file names. My best guess would be the input or output path, but I don't see the problem.


texts: How to retrieve Intersected geometry of the polygons geometry POINT (lat,lng) input?

texts:
I am passing a particular state (Indian states) latitude/longitude to the spatial query where it has to check intersection with the other table shapefiles polygon geometry. If that polygon is within the selected state polygons then it has to show all the state containing polygons.

Example situation I am trying is: suppose I have selected WEST BENGAL state its lat/lng will be passed and based on lat/lng intersection with other shapefile polygon has to done.

Working out on the spatial query

SELECT  distinct ST_AsGeoJSON(a.geom) As geometry ,b.proj_sno,b.title,b.start_date,b.end_date,b.total_project_cost,c.proj_category, a.state ,ST_Within(a.geom, ST_SetSRID(ST_GeomFromText('POINT(87.85 22.98)'),4326))  FROM proj_mst b,
 shape_file a,proj_cat_mst c,india_states d
 where  a.proj_sno = b.proj_sno  and
   a.proj_cat_id = c.proj_cat_id and a.proj_cat_id= 10

Structure for Understanding

1.

-- Table: shape_file

-- DROP TABLE shape_file;

CREATE TABLE shape_file
(
  gid serial NOT NULL,
  area_ha double precision,
  perimeter double precision,
  type double precision,
  objectid double precision,
  cname character varying,
  state character varying,
  geom geometry(MultiPolygon,4326),
  proj_sno integer NOT NULL,
  proj_cat_id integer,
  statecode integer,
  CONSTRAINT shape_file_pkey PRIMARY KEY (gid, proj_sno)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE shape_file
  OWNER TO postgres;

2.

-- Table: india_states

-- DROP TABLE india_states;

CREATE TABLE india_states
(
  gid serial NOT NULL,
  area double precision,
  perimeter double precision,
  indiasln_ integer,
  indiasln_i integer,
  state character varying(20),
  shape_leng numeric,
  shape_area numeric,
  geom geometry(MultiPolygon,4326),
  latitude double precision,
  longitude double precision,
  statecode integer,
  CONSTRAINT india_states_pkey PRIMARY KEY (gid)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE india_states
  OWNER TO postgres;

-- Index: india_states_geom_idx

-- DROP INDEX india_states_geom_idx;

CREATE INDEX india_states_geom_idx
  ON india_states
  USING gist
  (geom);

3.

-- Table: proj_cat_mst

-- DROP TABLE proj_cat_mst;

CREATE TABLE proj_cat_mst
(
  proj_cat_id serial NOT NULL,
  proj_category character varying(50) NOT NULL,
  CONSTRAINT proj_cat_mst_pkey PRIMARY KEY (proj_cat_id)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE proj_cat_mst
  OWNER TO postgres;

4.

-- Table: proj_mst

-- DROP TABLE proj_mst;

CREATE TABLE proj_mst
(
  proj_sno integer NOT NULL,
  group_id integer,
  div_id integer,
  proj_id integer,
  proj_cat_id integer,
  start_date timestamp(6) without time zone,
  end_date timestamp(6) without time zone,
  duration character varying(200),
  status_id integer,
  funding_agency_id integer,
  pi_ec_no integer,
  title character varying(200),
  title_abbrev character varying(20),
  study_area character varying(100),
  fa_role character varying(300),
  total_project_cost double precision,
  prog_sofar character varying(3000),
  prog_so_far_percent integer,
  project_type character varying(100),
  objectives character varying(3000),
  deliverables character varying(3000),
  tasks character varying(3000),
  nrsc_role character varying(300),
  proj_weekly_report integer,
  CONSTRAINT proj_mst_new_pkey PRIMARY KEY (proj_sno)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE proj_mst
  OWNER TO postgres;


texts: Layer not found in PostGIS db when it is listed by OgrInfo

texts:
Getting an error when trying to import a layer into R using the package RGDAL with the following:

geofile <- readOGR("PG:dbname=mydb", layer='temp_GIS_dw_file')

This returns the following error:

Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding,
use_iconv = use_iconv,  :
Layer not found`

Listing the layers with ogrListLayers shows the layer is present:

ogrListLayers("PG:dbname=mydb")

[1] "mb2013_clipped"    "mb_merged_testing" "temp_gis_dw_file"  attr(,"driver")
[1] "PostgreSQL"
attr(,"nlayers")
[1] 3

Furthermore, the following successfully imports the layer from the same DB:

geofile <- readOGR("PG:dbname=mydb", layer='mb2013_clipped')

I have cleaned up the geometry with ST_MakeValid() but this did not solve the problem.

What have I done wrong too late at night?

sessionInfo()
R version 3.3.0 (2016-05-03)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.1 LTS


texts: how to export the sublayers of a saved VRP layer to shapefile with arcpy

texts:
I have a bunch of saved layers of the result of VRP, I want to merge the route layers into one file, but apparently I first need to save them as shapefiles.
When I run this code for one of them, I get this error:

import arcpy
import os
import sys
from arcpy import env
directory = r'C:\\test'
shp_out = os.path.join(directory, 'output.shp')
layers = arcpy.mapping.ListLayers('results_MakeVRP.lyr',"Routes")[0]
arcpy.CopyFeatures_management(layers, shp_out)

Error:

Runtime error
Traceback (most recent call last):
  File "<string>", line 7, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\mapping.py", line 1517, in ListLayers
result = mixins.MapDocumentMixin(map_document_or_layer).listLayers(wildcard, data_frame)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 846, in listLayers
layers = self.layers
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 698, in layers
for frame in reversed(self.dataFrames):
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 710, in dataFrames
return map(convertArcObjectToPythonObject, self.pageLayout.dataFrames)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\arcobjects\mixins.py", line 694, in pageLayout
return convertArcObjectToPythonObject(self._mxd._arc_object.pageLayout)
AttributeError: 'str' object has no attribute '_arc_object'


texts: Identifying SBET coordinate system and reprojecting it with PDAL?

texts:
I've received discrete return LiDAR data from a data provider, which was acquired with an Optech C-300 scanner. In addition to the LAS data, the delivery includes the trajectory information in the Applanix SBET (Smoothed Best Estimate of Trajectory) binary format.

I was hoping to visualise the SBET data in QGIS - so I used PDAL SBET Reader to convert the file to ASCII as follows:

pdal translate input.sbet output.csv

The output file includes X, Y, Z coordinates, however these don't relate to my local projection (EPSG:2157):

"GpsTime" "Y" "X" "Z"
556975.005 0.933 -0.159 931.402
556975.010 0.933 -0.159 931.362
556975.015 0.933 -0.159 931.362
556975.020 0.933 -0.159 931.362
556975.025 0.933 -0.159 931.362
556975.030 0.933 -0.159 931.362
556975.035 0.933 -0.159 931.361

The min_x, min_y is (-0.1690, 0.9110) and the max_x, max_y is (-0.1650, 0.9130).

Using the SBET file, is it possible to identify which projection the X and Y coordinates relate to? And if so, can PDAL reproject these coordinates to EPSG:2157? Or is it something that needs to be corrected by the data provider?


texts: Displaying points from dataset within buffered areas of another using ArcGIS Online

texts:
I'm very new to using ArcGIS Online and I am creating a map for a project I have in class.  I have plotted all Chicago Public Schools and want to display JUST all arrests that occurred within a 5-minute walking distance of them, and hide the rest of the datapoints.  How can I do this?

I have tried creating drive-time areas and joining features but that did not work properly.

This is the link to my map: https://arcg.is/1mbr9j


texts:
I'm using qGIS 2.18.13 and I downloaded the plugin qGis2threejs to make 3D maps. I'm wondering how it works because it only works with some DEM rasters and not with everything. In particular, when showing the 3D map, it makes outputs like the one i posted.
I would also like to make 3D maps from Raster values not only with DEMs, but from other ones. Any suggestions?


texts: Snapping/moving points near polygon to inside polygon using QGIS

texts:
I have a set of polygons (buildings) and points that I'd like to spatially join. Unfortunately, some of these points are positioned outside the polygon because all points are positioned on the centroid of each building and sometimes when a building is shaped like an L, for example, the point feature isn't placed inside the polygon.

Now I'd like to preferably use QGIS to move the points that are outside polygons to be situated inside the nearest polygon.

I've tried the Snap Geometries to Layer function but the points sometimes end up on a node which is a corner node of more than one building.

Any ideas?

Grey buildings are the ones where the point is positioned outside the building boundaries and therefore can't be joined.


texts: GRASS calculation of centroids : very important differences with QGIS calculation

texts:
I use v.extract module in a python script :
    grass.run_command('v.extract', input=UH, output=centroids, type='centroid',
                                  overwrite='True', stdout = DEVNULL, stderr =      DEVNULL)

when I compare the result with result from calculation with QGIS I notice some differences.

red : GRASS v.extract, blue : QGIS.
One of my colleague observed the same with its own shapefile.
SAGA module (Polygon centroids) gives same result as QGIS.
It seems that v.extract calculation isn't correct.


texts: Extracting polygons from road network using ArcPy?

texts:
I need to extract faces (or polygons) from my road network and calculate their areas.

I mainly use ArcPy.

If it needs to be converted to graph structure, include it as well.


texts: What is the unit of "output cell size" in IDW in ArcGIS?

texts:
I have a wind speed point data which I interpolated using IDW. While interpolating, it asks for output cell size. I am confused as in it is in which unit?

My data is in UTM projection and WGS 84 datum. Map units and display units are in meters.

I am using ArcGIS 10.2.2.

Reason that I am confused because if I want 5 Kilometer cell size and I give 5000 value to "output cell size", then it is not taking the value and error is coming.


texts: What is the most convenient way to delete geometries from a rapid eye image? I want to draw a polygon and delete that part from the image. Like the opposite to a clip function, that keeps the geometry but not the rest of the image.


texts:
I'm having trouble exporting a map from QGIS to an image on Linux. The exported image has much stronger color contrast than inside QGIS. This occurs both on PNG and TIFF exports.

I exported by Project > Export/Import > Export Map To Image.

On the left is the exported TIFF (PNG looks the same). On the right is the map in QGIS. Notice how much darker some of the green areas are on the left image.

Is this a known issue? Is there a better way to export an image from QGIS?

EDIT: After playing around with this a while longer, it seems that there is a relationship between the size of the exported image and the degree of constrast/darkening. If I export with only 3k pixels to a side, the image looks like it does in QGIS. As I increase this, it starts to change more and more.


texts: Intersection between polylines feature to polygon in QGIS

texts:
I would like to transform a line shapefile into polygons. However, the tool Vector>Geometry>Line to Polygon doesn'to work to me. What I would like to do is using the intersection between lines and make them as vertices of a polygon.
In ArcGis, the tool Feature to Polygon does exactly what I want.
Does someone know a way of doing so in QGIS?


texts: Calculating Surface Roughness from Land Cover Data/Land Elevation

texts:
I'm trying to create a surface roughness file over the continental United States in NetCDF format. I've surface elevation file of 200 meters resolution and land cover data of 30meters resolution. I'm trying to create the roughness file in ArcMap.

Can anyone suggest any approach which can be used to calculate the roughness?


texts: I have two feature classes, one ploylines, one a buffer.  How can I populate a field of the polylines to indicate if they intersect?


texts: Reprojecting data with GeoPandas using 'to_crs'?

texts:
I am trying to reproject a polygon so that it has the same projection than a raster dataset.

raster.crs
CRS({'proj': 'tmerc', 'lat_0': 0, 'lon_0': -183, 'k': 0.9996, 'x_0': 500000, 'y_0': 0, 'datum': 'WGS84', 'units': 'm', 'no_defs': True})

polygon.crs
{'init': 'epsg:3035', 'no_defs': True}

I believe I should be using an instruction like the one below:

polygon = polygon.to_crs(epsg=XXXX)

where XXXX is the epsg code that tells the coordinate system of the raster dataset. However, I don't know hot to obtain such code. (I know that the epsg code of WGS84 is 4326, but I would like to know how to obtain this number programmatically)


texts: Snapping vectors to vertices for ROI delineation?

texts:
I am trying to manually extract a set of ROIs in ENVI 5.0.

However, I can't seem to successfully line up the vector to be flush with the boundary of the pixels, which means I don't have fine-tune control over which pixels are included/excluded at the edge of my vector. Also, the resulting ROI ends up splicing pixels in half.

Is there anyway I can snap the vertices to the pixels when I'm drawing my vectors/ROIs?


texts:
I have a table made of merged point datasets in PostGIS.  I want to select all of the features where two or more points have identical geom values.

At present I just want to select so I can review the associated field values.


texts: How to count the number of data points in each area in QGIS?

texts:
I have a shapefile (.shp file) which defines a large number (several thousand) separate areas.

I also have an Excel spreadsheet which contains a large number of data points. Each data point consists of a latitude and longitude, and a "weight".

I want to determine, for each area in the shape file, the total weight of all the data points whose lat/long is in that area. I would prefer to have this information in a format that can be downloaded back into the Excel spreadsheet.

What is the best way to do this in QGIS?


texts: Shapefile/vector layer not appearing in QGIS

texts:
I have a shapefile (calling it 'new.shp' for discussion's sake) that won't show up with other layers.

If a fresh file is started and I bring in new.shp, it will appear, but any other layers imported after does not. On the flip side, if other layers are imported first, the other layers appear, but new.shp will not show up.

That's not too much detail to troubleshoot from, but I was hoping someone might be able to suggest what the issue might be?

I've just started learning GIS and using QGIS, so I'm a bit lost at the moment.

Other notes:
- I looked into the CRS for new.shp which says (EPSG:4326, WGS 84); The other layers are also set to this.
- Using QGIS 2.10.1 on Windows 7


texts: Editing attributes in ArcMap causes conversion of arcs to densified lines unintentionally

texts:
I am using postgresql GDB for storing and managing polyline feature classes through arcmap but when I create arc feature along with filling any attributes the arc is being converted to a densified straight polyline unintentionally.

It is noticeable that if I do not fill attributes when editor is started, arc structure (just 2 end points) is being preserved and saved as well.

I found that gdb_geomattr_data column is a repository for curve parameters but my feature class (which has been registered as versioned) not showing this trend! I got completely confused.


texts: How to use dateRangeContains in Google Earth Engine?

texts:
I need the closest image of a collection to a date. Let's say the image that is closest to 1st of August 2019.
The help says that
Creates a unary or binary filter that passes if the left operand, a date range, contains the right operand, a date.

which seems to be what I need.
I expect to use it like this
ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filter(
    ee.Filter.dateRangeContains(ee.dateRange('2019-07-24','2019-08-09'),'2019-08-01')
  )

but the problem is that dateRangeContains has the following arguments:
ee.Filter.dateRangeContains(leftField, rightValue, rightField, leftValue)

leftField (String, default: null):

A selector for the left operand. Should not be specified if leftValue is specified.
rightValue (Object, default: null):

The value of the right operand. Should not be specified if rightField is specified.
rightField (String, default: null):

A selector for the right operand. Should not be specified if rightValue is specified.
leftValue (Object, default: null):

The value of the left operand. Should not be specified if leftField is specified.

Considering the definition, I did left the first 2 arguments aside, and used only the last 2, which are named in the definition, i.e. the left operand, a date range, contains the right operand, a date:
ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
  .filter(
    ee.Filter.dateRangeContains(undefined,undefined,ee.dateRange('2019-07-24','2019-08-09'),'2019-08-01')
  )

Needless to say, it didn't work, although there is an image with the date on 2019-07-25


texts:
I want to extrude a polygon shape in postgis to create a pseudo 3D effect. To this end I have written a crude function to achieve it. This is very much test code and it creates a new Y vertice for each point on the polygon and then closes it by returning to the original point:-

CREATE OR REPLACE FUNCTION public.extrude_polygon(wkb_geometry_param geometry, height integer, simplify boolean DEFAULT false)
  RETURNS geometry AS
 $BODY$
DECLARE
f int;
ret_geom geometry;
wkb_geometry geometry;
BEGIN

--convert polygon to linestring
IF ST_GeometryType(wkb_geometry_param) != 'ST_Polygon' THEN
    RETURN NULL;
END IF;

IF simplify THEN
    wkb_geometry =         ST_Simplify(ST_Transform(ST_Exteriorring(wkb_geometry_param), 27700), 0.5);
ELSE
wkb_geometry = ST_Transform(ST_Exteriorring(wkb_geometry_param), 27700);
END IF;

--initialise output geometry
ret_geom =ST_MakeLine(ST_PointN(wkb_geometry,1),ST_PointN(wkb_geometry,1));

--Move first point to up

SELECT ST_AddPoint(ret_geom,
      ST_MakePoint(ST_X(ST_PointN(wkb_geometry, 1)),
                   ST_Y(ST_PointN(wkb_geometry, 1)) + height)
) into ret_geom;

FOR f IN 1..ST_NPoints(wkb_geometry) LOOP
    IF f < ST_NPoints(wkb_geometry) THEN
    --across to next high point
    SELECT ST_AddPoint(ret_geom,
               ST_MakePoint(ST_X(ST_PointN(wkb_geometry, f + 1)),
               ST_Y(ST_PointN(wkb_geometry, f + 1)) + height)
    ) into ret_geom;
    --down to next point
    SELECT ST_AddPoint(ret_geom, ST_PointN(wkb_geometry,f + 1)) into ret_geom;
    --back to last point
    SELECT  ST_AddPoint(ret_geom, ST_PointN(wkb_geometry,f)) into ret_geom;
    --back then up again
    SELECT ST_AddPoint(ret_geom, ST_PointN(wkb_geometry,f + 1)) into ret_geom;
    SELECT ST_AddPoint(ret_geom,
               ST_MakePoint(ST_X(ST_PointN(wkb_geometry, f + 1)),
               ST_Y(ST_PointN(wkb_geometry, f + 1)) + height)
    ) into ret_geom;
ELSE
    --across to first high point
    SELECT ST_AddPoint(ret_geom,
               ST_MakePoint(ST_X(ST_PointN(wkb_geometry, 1)),
               ST_Y(ST_PointN(wkb_geometry, 1)) + height)
           ) into ret_geom;

    SELECT ST_AddPoint(ret_geom, ST_PointN(wkb_geometry,1)) into ret_geom;

END IF;
END LOOP;

RETURN ST_Buffer(ST_Buffer(ST_MakePolygon(ret_geom),10), -10);

END;

$BODY$
  LANGUAGE plpgsql

It works with simple polygons but has problems with interior rings, but the main issue is that it is really slow. I need to output the resultant shape as a polygon that can be shaded and rendered in mapserver. Hence the buffer operations at the end which is the only way I know of reducing the shape to it's outline.

The end result will be an extruded shape representing the original polygon. I can then offset the original polygon by the same extrusion distance and place it on top to make the roof.

I considered using the ST_Extrude function in postgis-2.1.1 BUT this creates a ST_PolyhedralSurface type and I am not able to render it in mapserver. As far as I can tell there is no way to create an outline of this either as ST_Buffer does not work with ST_polyhedralsurfaces.

So, my question is, can my function be improved? Or is there a better approach. The output needs to look as per the diagram which I created by placing the offset polygon onto my extruded shape.


texts: Shapefile Layer display filtering in QGIS?

texts:
I use QGis 3.0, created several layers, each with multiple attributes.

I want to display only objects which have f.e. "123456789" as value in a specific attribute to print it and repeat with the next value.
It shouldn't just select the items but make everything else invisible.
I don't want to filter every singly layer. A filter which applies to every layer would be perfect since the attribute has the same name and value in every layer for the same order number.

The list might make thing a bit more clear:

The upper layer group "Trasseninhalte" are cables and pipes and connections.

The lower side "Trassenobjeckte und APL's" shows burrows to be made or already existing, usable pipes. "Knotenpunkte", "APL" are just connection points.


texts: Changing ECW output colour for no input value areas

texts:
I am attempting to convert a mosaic geodatabase raster to ECW using FME, but get black image output around the edge where there is no input value.
How do I change this from black?


texts: What to do with unused layers in my personal geodatabase

texts:
I have a question about my personal geodatabase. You can see it at the picture below.

I started with the database in Arcmap with my normal data (Black sea, Danube, East Aegean and West Aegean). In access I added the Excel table with data and linked that data with the layers by the use of queries.

So in Arcmap I'm using the query_.... layers, but the normal ones are still there. Does somebody know if there is a way or if it is possible to remove  the old layers, without this will influence the Query which is half based on the old layers?

I already tried to delete the Records of the old layers in the GBD_GeomColumns table in mc Access, but they still stay visible at the personal geodatabase in arc catalog (I refreshed it). When I delete for example Black Sea the query? Black_sea isn't working anymore because mc Access can't find the normal Black Sea table.

So maybe there is a way how the old layers stay at the database, but aren't visible at ArcCatalog?

Hope someone knows how it works, because it would be a lot more synoptic in use.


texts: QGIS and batch processing: Is there a procedure to fill multiple parameters using arrays?

texts:
I'm trying to calculate hillshades (Raster analysis -> Hillshade) for different time of the year, using batch processing and different input values for Azimuth and Altitude of light.
I figured out how to set the parameter Azimuth of the light using “Add value by expression...“ [dummy exemple: array(45, 50, 55, 60)]. I get four new lines filled with the expected azimuth values and the default value for Altitude of the light (45.0).
If I then set the altitude following the same procedure (“Add value by expression...“, array(100, 110, 120, 130)), I get 4 new lines filled with expected values for Altitude and default value for Azimuth (315.0).
Is there a procedure to set both Azimuth and Altitude in each line using arrays?


texts: How to add start and end coordinate to attributes of a road network?

texts:
I have generated a road network from open street map data. How can i add start and end coordinates (point_x, point_y) to its attributes table via toolbox or arcpy?

That's how the roads look like.
I also have the shapefile, but i believe this should be more appropriate, since the junctions i need are already identified.


texts:
I'd like to show the user the current map scale. I update the current map scale in a textblock  everytime the map size/scale has changed. This is my code, to show the map scale (in C#):

private void Map_SizeChanged(object sender, SizeChangedEventArgs e){
    MapScaleTextBlock.Text = String.Format("Scale {0:N}", Map.Scale); //Note: It has to be formatted
}

and here is the xaml code:

<TextBlock Name="MapScaleTextBlock"></TextBlock>

The output is always this:

  Scale n. def.

Whats wrong? I searched already but I found nothing helpful.

I'm using ESRI API for WPF and C# in Visual Studio 2010.


texts:
I have been trying to mask clouds with the quality control band with landsat 8 imagery.  The problem I continue to run into is that the quality control bands mistake tops of buildings and concrete as clouds.

Have you guys had better success with Fmask?  I tried using it earlier, but couldn't find out where it put the file.  The website said it would be called XXXFmask.  Would I even be able to open Fmask files in Arcmap?

These are the steps I've taken.
1) Built raster attribute table for QC band (see below).
2) Selected only the pixels whose values said there were clouds present.
3) Used Raster to Polygon tool to convert selected pixels to a shapefile.
4) Selected all the cloud polygons and used mask tool on image analysis window to remove data from the cloud cells.
5) Now I plan to run my image classification tools.  I've included an image of Greenville, TN before and after the raster is masked.  NoData cells are displayed in pink for visibility.


texts: PostGIS ST_Area with similar geographies producing different results

texts:
I'm running ST_Area on two different NaturalEarth datasets but am getting different results when the geographies are clipped before performing the area calculation.

First, I've loaded the 10m "admin0" dataset (http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip) into a table named countries.

Second, I've loaded the 10m "land" dataset (http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/physical/ne_10m_land.zip) into a table named land.

The area calculation is basically the same for both:

SELECT sum(ST_Area(geom::geography)) AS area FROM land;
     area
-----------------
 146826440837593

SELECT sum(ST_Area(geom::geography)) AS area FROM countries;
       area
-----------------
 146826603523460

Examining the data, it appears that both datasets are the same except that the land table has one large multipolygon, whereas countries has many individual multipolygons (ie one for each country in the world).

The issue arises when I clip both to a latitude of +/-52 degrees, in which case the areas I calculate for each diverge quite a bit. This is the query:

SELECT
sum(
  ST_Area(
    ST_Intersection(
      ST_GeomFromEWKT(
        'SRID=4326;POLYGON((-180 -52, -180 52, 180 52, 180 -52, -180 -52))'),
      geom
    )::geography
  )
) AS land_area
FROM
 land;

And it produces these results:

    land_area
-----------------
 119910798321370

  countries_area
-----------------
 106619873643318

Any ideas as to why this would be the case? I've also visualised the clipped geometries for both tables and they appear the same to me and I don't see the same area differences when I clip to latitudes like +/-80 degrees.

All the geometries appear to be valid before the intersection as well as after the intersection. For example, the following query produces no results:

SELECT
 *
FROM
 countries
WHERE
NOT ST_IsValid(
      ST_Intersection(
        ST_GeomFromEWKT(
          'SRID=4326;POLYGON((-180 -52, -180 52, 180 52, 180 -52, -180 -52))'),
        geom
      )
    );

UPDATE: If I do the calculation without converting to geography, the numbers do not diverge like they do above:

    land_area
------------------
 9942.37708696352

   country_area
------------------
 9942.38832452307


texts: How to escape percent in a LIKE expression in qgis?

texts:
How can we escape a percent character in a LIKE expression in QGIS?

For example, if I have the following column:

123% , 15% , albert, greg, 15%

How do I select all the values with % in it without changing the field?

It seems possible to do this in ArcGIS, according to this answer but I have not found an equivalent solution with QGIS.

I have tried escaping with %%, \%, \\\%, [%], ^%, etc. but without success.

[EDIT]

I actually truly was wondering if it is possible to escape the % (and _) character, rather than to use perfectly (valid) work-around. The example provided was but for illustration.
Any ideas?


texts: Getting latLngBounds for tile based on its tile coordinate using Leaflet

texts:
I would like to download raster image data in a leaflet / esri-leaflet context.  To provide a little backround info: When you add an esri-leaflet raster layer to a leaflet map, you can see that on map move end, a request is made to the arcgis layer service url.  The request contains the bounds of the map, and an image of just those bounds are returned and layed over the map as an image overlay.
I would like to tweak this process. (I am reproducing this behavior on a node server, without a map.) Instead of requesting one single image that covers the whole map bounds, I want to download a series of tiles, which when put together, make up the same image that covers the whole bounds of the map.  I would like these tiles to have the same bounds as as the tiles of a standard slippy map tile layer (assuming standard web mercator projection).
Let's say I have a function already written that takes in a latLngBounds, and returned the fetched image for those bounds.  (Its too complicated to post it all here, as it is reconstructed from a lot of esri-leaflet's source code).  We'll call it fetchImageFromBounds(bounds: L.LatLngBounds).
For a given map latLngBounds, I can get all the tiles contained within those bounds using xyz-affair:
import xyz from 'xyz-affair';

const scale = 10;

const fullExtentBounds = [topLeft, bottomRight];
const tileCoords = xyz(fullExtentBounds, scale);

Now I would like to write a function that returns the latLngBounds for each of the tile coordinates:
function tileCoordToBounds(coord){
  const {x, y} = coord;

  const topLeftLayerPoint = new L.Point(
    X * tileSize,
    Y * tileSize
  );

  const bottomRightLayerPoint = new L.Point(
    (X + 1) * tileSize,
    (Y + 1) * tileSize
  );

  return L.latLngBounds(
    L.CRS.EPSG3857.pointToLatLng(topLeftLayerPoint, scale),
    L.CRS.EPSG3857.pointToLatLng(bottomRightLayerPoint, scale)
  );
}

const tilesBounds = tileCoords.forEach( coord => {
  const singleTileBounds = tileCoordToBounds(coord)
  fetchImageFromBounds(singleTileBounds)
});

I recognize that the bounds of tiles contained within fullExtentBounds will reach outside of fullExtentBounds.  I'll end up with a larger bounds than I started, and that is fine.  Here's a visualization:

I want to make sure that these bounds produce a series of images that is continuous, meaning there are no gaps between tiles, and no overlaps, down to the exact pixel of the original raster source.  Meaning that I could lay these tiles out like in the diagram, and they would produce exactly the same image as fetchImageFromBounds(fullExtentBounds).
Does tileCoordToBounds look right, in that the array of bounds are continuous, meaning without gaps and without overlaps?


texts: Leaflet.locate: permission error on different browsers

texts:
I inserted a simple button for geolocation into a Leaflet-generated map by adding

L.control.locate().addTo(map);

into my HTML file.
In Chrome and Firefox it works great, Safari and Internet Expolorer are rejecting the function by showing following error:

  english: "geolocation error: This website does not have permission to use the Geolocation API"
  german: "geolocation error: Diese Website verfügt nicht über die Berechtigung zur Verwendung der Geolocation-API"

Do you have any idea, where this issue can come from?
I've tried to use the functionality on different browsers on different devices (desktop and mobile), all with the same error. I've checked the location permissions on the devices as well.


texts: Persisting ArcGIS Network Analysis Route Layer?

texts:
I'm using ArcGIS 10.0.
I'd like to be able to specify the location (path) where to create the route layer, for network analysis, when I click on Network Analysis -> New Route. Is that possible?

If it is not, how can I have access to the location where the route layer is created?

I'm using the route layer in a task of a geoprocessing service, and I'd like to give to the task the exact location (absolute path) of the route layer.

Could anybody help?


texts: QGIS fails to recognize Geometry Column in SQL Server

texts:
I have a table [stg].[Stops] in SQL Server that stores a variety of fields including a well-known binary, the GEOM column, that has been created from the stop_lat and stop_lon fields. The GEOM is a geometry data type. Here is a screenshot of the table.

However, when I go into QGIS and attempt to add the SQL Server table. QGIS does not recognize the table in the [stg] Schema.

What am I missing here? I do not have a geometry_columns table, but I don't think thats the issue.


texts: AttributeError when using arcpy.GetParameterAsText()?

texts:
I found a script to batch export geotiffs from data view and when hard coded it work well

mxd = arcpy.mapping.MapDocument("CURRENT")

res = 200 #resolution of tif
export = "S:\\6 - GIS\\For Randoms\\For Jeff\\Thunderchild Reserve\\Shapefiles\\Imagery" #folder
preFix = "image" #filename
widthcolumn = 3910 #width (number of columns)
heightrow = 1765 #height (number of rows)

df = arcpy.mapping.ListDataFrames(mxd)[0]

for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
  mxd.dataDrivenPages.currentPageID = pageNum
  arcpy.AddMessage("Exporting Map " + str(pageNum) + " of " + str(mxd.dataDrivenPages.pageCount))
  arcpy.mapping.ExportToTIFF(mxd, export + "\\" + preFix + str(pageNum) + ".tif",df,
                 df_export_width=widthcolumn,
                 df_export_height=heightrow,
                 resolution = res,
                 world_file=True,
                 tiff_compression="NONE",
                 geoTIFF_tags=True)

When I try to define my variables using the arcpy.GetParameterAsText() and set my export variable as folder in the tool box I return an error.

mxd = arcpy.mapping.MapDocument("CURRENT")

res = arcpy.GetParameterAsText(0) #resolution of tif
export = arcpy.GetParameterAsText(1) #folder save path
preFix = arcpy.GetParameterAsText(2) #filename
widthcolumn = arcpy.GetParameterAsText(3) #width (number of columns)
heightrow = arcpy.GetParameterAsText(4) #height (number of rows)

df = arcpy.mapping.ListDataFrames(mxd)[0]

for pageNum in range(1, mxd.dataDrivenPages.pageCount + 1):
  mxd.dataDrivenPages.currentPageID = pageNum
  arcpy.AddMessage("Resolution = " + res)
  arcpy.AddMessage("Export Folder = " + export)
  arcpy.AddMessage("Exporting Map " + str(pageNum) + " of " +     str(mxd.dataDrivenPages.pageCount))
  arcpy.mapping.ExportToTIFF(mxd, export + "\\" + preFix + str(pageNum) + ".tif",df,
                 df_export_width=widthcolumn,
                 df_export_height=heightrow,
                 resolution = res,
                 world_file=True,
                 tiff_compression="NONE",
                 geoTIFF_tags=True)

The error message that I get is:

*Executing: BatchExportGeotiffs 300 "S:\6 - GIS\For Randoms\For Jeff\Thunderchild Reserve\Shapefiles\Imagery" image 4888 2647
Start Time: Fri Mar 03 16:19:40 2017
Running script BatchExportGeotiffs...
Resolution = 300
Export Folder = S:\6 - GIS\For Randoms\For Jeff\Thunderchild Reserve\Shapefiles\Imagery
Exporting Map 1 of 4
Failed script BatchExportGeotiffs...
Traceback (most recent call last):
  File "S:\6 - GIS\Resources\scripts\Python\BatchExportGeoTiffs.py", line 22, in <module>
    geoTIFF_tags=True)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1391, in ExportToTIFF
    layout.exportToTIFF(*args)
AttributeError: DataFrameObject: Error in parsing arguments for ExportToTIFF
Failed to execute (BatchExportGeotiffs).
Failed at Fri Mar 03 16:19:40 2017 (Elapsed Time: 0.04 seconds)*

Can anyone tell me where I'm going wrong?


texts: Converting TIR data (CSV) to Raster layer using QGIS

texts: I have thermal infrared data from a FLIR camera. I would like to be able to import this to Q with spatial reference and temperature data as a raster. Does anyone have experience with flir radiometric images?


texts:
Consider the following code:
import numpy as np
from shapely.geometry.polygon import LinearRing

llpts = np.asarray([[18.5165    , 73.82883333],
       [18.51678333, 73.8279    ],
       [18.5156    , 73.82848333],
       [18.51586667, 73.82755   ],
       [18.5165    , 73.82883333]])

linring_ll = LinearRing(llpts)

print('linring_ll.is_valid =',linring_ll.is_valid)

Using Shapely 1.7.1 and Python 3.9.2 the above code prints 'True' on my computer even though it seems that the LinearRing crosses itself (I plotted it in Spyder). Is this a rounding / numerics issue? This goes away when I convert the points from lat/lon space to xy space using equirectangular projection.
The complete code follows:
import numpy as np
from shapely.geometry.polygon import LinearRing, LineString, Polygon

def convert_llpts_xypts(latlon,rad_earth,center_lat,projtype='equirectangular'):
    # convert from latlon space to xy space
    # latlon must be in degrees
    # lat lon looks like [(lat1,lon1),(lat2,lon2)....]
    # lat lon needs to be a np array

    if (projtype == 'equirectangular'):
        # https://stackoverflow.com/questions/16266809/convert-from-latitude-longitude-to-x-y
        _radconv          = np.pi/180.0
        _rade_radconv     = rad_earth*_radconv
        _rade_radconv_cos = _rade_radconv*np.cos(center_lat*_radconv)

        _xx = _rade_radconv_cos*latlon[:,1]
        _yy = _rade_radconv*latlon[:,0]
        _xx = _xx[:,None]
        _yy = _yy[:,None]

        return np.hstack((_xx,_yy))

    else:
        raise ValueError(f'Invalid projtype={projtype} in convert_llpts_xypts')

# points in lat/lon space
llpts = np.asarray([[18.5165    , 73.82883333],
       [18.51678333, 73.8279    ],
       [18.5156    , 73.82848333],
       [18.51586667, 73.82755   ],
       [18.5165    , 73.82883333]])

# compute a latitude close to the center of the polygon
center_lat = np.sum(llpts[0:4,0])/4.0

# convert points to xy-space
xypts = convert_llpts_xypts(llpts,6378137.0,center_lat)

# make LinearRings
linring_ll = LinearRing(llpts)
linring_xy = LinearRing(xypts)

print(f'validity of ring in lat/lon space = ', linring_ll.is_valid)
print(f'validity of ring in xy space = ', linring_xy.is_valid)


texts:
I am real new to C++ and am having trouble compiling a C++ program that uses the proj.4 library. I downloaded and compiled and installed the proj.4 library on my machine. So, the /usr/local/lib'contains

libproj.a
libproj.la
libproj.so -> libproj.so.0.6.6
libproj.so.0 -> libproj.so.0.6.6
libproj.so.0.6.6

and the /usr/local/include contains

nad_list.h
org_proj4_Projections.h
proj_api.h
projects.h

and to compile my program test.c, I am using the following command

g++ -Wall -g first.c -lproj -o test

The compilation goes fine. But when I try to execute the program with ./test I get the error message

./test: error while loading shared libraries: libproj.so.0: cannot open shared object file: No such file or directory

I can't tell what I am missing here. Any help or suggestions here would be very helpful.


texts: How to aggregate points with same value into polygons from a shapefile using GDAL or any other opensource solution

texts: I have a shapefile with around 19,000 points. Its basically export from a raster. No i need to extract polygons, by aggregating the points which have same value.The field who's value i am going to use for aggregation is dynamically calculated each time using the elevation of points. NOw i need to spit out polygons. How can I do that using GDAL? is there a utility to do it. Any other opensource solutions are welcome. I have ArcGIS which has a toolbox called 'Aggregate Points' but somehow licence for it is missing.


texts: replacement for cartodb.SQL({user: 'user'}).execute(sqlquery) in new cartojs

texts:
this used to be my code in older carto.js

cartodb.SQL({user: 'user'}).execute(sqlquery).done(function (data){
  console.log(data);
}

how to readapt it in newer cartojs 4.0 beta


texts: 3D modelling mine plans - how to set elevations?

texts:
I am trying to make a 3D model from an old mine plan, I have created a new project in QGIS and added a new geo package layer and checked include Z dimension. Then I have traced a line from the old plan with the add line feature. In vertex editor I have entered the Z values as -10 for all points.
But when I look in the 3D viewer all the lines still seem to be on the zero plane
Am I missing something like Z is not elevation? Or am I going completely wrong with this?


texts: DEM disappears at certain zoom level in browser

texts:
my Qgis2threejs-model looks fine in the plug-in window, but in the browser (e.g. Edge) I have to zoom out to see the DEM. If I zoom in the DEM desappears at a certain zoom level. I don't want this behaviour.
Is there a parameter in any of the js-files, to define this zoom level or even to switch it off?

QGIS 3.4.13 and qgis2threejs 2.3

Kind regards!

]2


texts: Setting Max Cache for GDAL, using C# bindings

texts:
The C# bindings in GDAL expose the method called

  Gdal.SetCacheMax(int nBytes)

for setting the GDAL Raster block cache size. However, it takes in an int and that int represents the number of bytes. So, to set the memory to, let's say, 8GB I will have to provide 8 * 1024 * 1024 as the input value, which exceeds the maximum allowable int size and wraps around as a negative number.

How do I set like 75% of my RAM OR like 12 GB of RAM for CACHE_MAX value using C# bindings of GDAL?


texts: Using complex overlapping polygon in QGIS?

texts:
I'm trying to solve quite a tough analytical problem and hoping I can find a way to do it in QGIS. Basically, I have developed a map of Australia with an overlay of agricultural production by zone (polygons linked to a dataset of tonnes of grain produced per zone). In addition, I have mapped all of Australia's railway lines in a vector layer. What I need to do is calculate what fraction of the grain production is within X km of a railway line.

The mathematical formulation means I effectively need to draw a polygon of radius X km around all my railway lines and then calculate the area of overlap with the agricultural polygons multiplied by the production intensity of those overlap zones.

Is it at all possible to do something like this with QGIS?


texts: How to set up an ArcMap Relate with a table from MS SQL Server 2012

texts:
This is my first time working with data directly from MS SQL Server 2012, and I'm unsure whether I can use it in a Relate in ArcMap.

I have a table in MS SQL Server, which I brought into ArcMap via the Database Connections option. When I add this table to ArcMap using the standard Add Data button, a New Query Layer dialog appears:

ArcMap requires a unique identifier field, however this table doesn't have one. I chose one of the existing fields, even though it's not unique (the first problem?).

I can now view the table in ArcMap, and I'm able to set up a Relate between this table and a featureclass, based on the matching BoreholeID field.

However, when I select a borehole from the featureclass, or a row in the table, the Relate is displayed but greyed out:

Can I use a table from MS SQL Server Management Studio in an ArcMap
relate?
If so, is the cause of the problem the lack of a unique ID
field?
If not, what else could be causing the problem?


texts:
I was wondering if its possible to make clips and other types of operations with a shp or tiff file, but all worked on a csv as using text?
What I want to do is replicate these operations in labview in csv.

I am working on a project, with shape files and raster files, and I was looking at a website that converts csv, among other formats to shape and raster files. What I want to know is the math involved when doing a clip for example, so I can replicate that in a csv file inside labview. I just want to know the math behind the processes in qgis or arcgis, so I can replicate these processes in labview, with a csv file.


texts: How to insert a quotation mark through the Field Calculator of attribute table of a Shapefile?

texts: I have a shapefile of polylines with information of diameter. But the problem is that i didn't enter the inch sign(") with the values of diameter. What would be the solution that qutation mark would place along with diameter values?


texts: Memory leak updating geometry - ArcGis Runtime .Net

texts:
i'm working with ArcGis Runtime SDK for .Net v10.2.5.

I have an UDP socket listening and waiting for image data that fires a function executed on a different thread in background.

I want to draw the image over a ellipse of arbitrary radius so i use

var filestream = System.IO.File.Open(imagepath, FileMode.Open, FileAccess.Read);
MapPoint point = new MapPoint(center.longitude, center.latitude, SpatialReferences.Wgs84);

var polySymbol = new Esri.ArcGISRuntime.Symbology.PictureFillSymbol();
await polySymbol.SetSourceAsync(filestream);

var param = new GeodesicEllipseParameters(point, 25, LinearUnits.Meters);
var ellipse = GeometryEngine.GeodesicEllipse(param);

***//HERE IS THE PROBLEM***
_graphicsLayer.Graphics.Clear();
_graphicsLayer.Graphics.Add(new Graphic { Geometry = ellipse, Symbol = polySymbol });

This is done ~5 times/second. Despite i'm clearing the layer each iteration there is a memory leak that is increasing memory use till app crashes.

I read about problems with memory using ArcGIS and Geometry process, so i'm not sure if i'm hitting a wall or just doing things badly


texts: Saving Label Text Size, Placement, Wrapping, As Template using QGIS

texts:
I am using QGIS 3.4.

Is it possible to export & save formatted labels as a template in a file that can be opened in future projects?

Ideally, it would save my manual label placement, text size, wrapped characters/alignment to save time in future projects.  I figured out how to do this for a customized color ramp, but cannot figure out how to do it for labels.


texts: Icon inside Leaflet control layer does not display?

texts:
I created a leaflet control layer to switch between two base maps (Bing and OpenCycleMap), Yet I don't see the icon inside the layer control:

var baseMaps = {
             "Bing": bingLayer,
              "OpenCycleMap": tileLayer
         };
        L.control.layers(baseMaps).addTo(map);

The result in the navigator:

Is there any special config to add the layer icon instead of the blank space? As I know this icon should be by default displayed.


texts: Visualizing GeoServer web map in other computers of company?

texts:
I made a Leaflet web map that I have no problem visualizing in my computer. It is a simple web map that contains a basemap and a few shapefiles (it's only use is to visualize the spatial data and when clicking, show some attributes), served from Geoserver on my computer.

Everything works just as intended, but now I want other people in my company to be able to use it. The problem is that when they enter the web map, they see the basemap ('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png') but not the shapefiles, so it is clear that the shapefiles are not in the right place.

I am using Geoserver with Tomcat, and I have the shapefiles in my computer. The data directory in Geoserver is "C:\Program Files (x86)\GeoServer 2.14.1\data_dir", where the shapefiles are.

My guess is that when they enter the web map, they don't have access to my "C:\Program Files (x86)\GeoServer 2.14.1\data_dir" and that is why they don't see the shapefiles.

Is there some way to save the shapefiles in geoserver/tomcat so everyone can visualize them, or another way to deal with this?

It sounds like some basic knowledge, but I am a beginner with geoserver, tomcat and webmapping and I have been searching for days without success.

I have also tried to put the shapefiles in a shared folder and set the shapefile directory of the Geoserver Store as that folder, but that doesn't solve the problem either.

Next, I'm showing part of the code, as required in the comments. Everyone sees the OpenStreetMaps basemap, but no one except me sees the 'Mun' layer, from Geoserver.

var osm = L.tileLayer('http://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    attribution: 'OpenStreetMaps',
    maxZoom: 18
}).addTo(map);

var mun = L.tileLayer.betterWms("http://" + serv + ":8089/geoserver/mex/wms?", {
    layers: 'Mun',
    format: 'image/png',
    transparent: true,
    opacity:0.3,
    version: '1.1.0',
    attribution: "Geoserver"
    }).addTo(map)


texts: Show attributes of area of interest by using Atlas

texts:
I use the "Atlas" tool of QGIS 2.18 to generate maps of different districts of my "city" layer. The information shown inside every district are the name, a value for power, and a pie chart. When I change between the single maps of the atlas, I would like to have only the information shown, which are interesting for the special district. The information of the residual district should be hidden.

I already found a way to select only the district of interest by using $id = $atlasfeature for stile and label, but I can't hide the  diagram...

Is there a way to do so?


texts: Routing for Waste Collection (U-Turn problem)

texts:
A general question about routing problem: I am working on waste collection problems and I am having some trouble with some routing problems (the truck has to pass by every bin). So far I have tried a few things:

•   GRASS traveling Salesman solver (in QGIS and GRASS directly)

•   GVSig traveling salesman but not with better results.

•   GrassHopper (https://www.graphhopper.com/)

TSP in GRASS seems to work fine is you have few points and relatively large distances between points. But because I have a large (100-300 points in close proximity) the path created is not realistic, mostly due to U-Turn problems (the shortest path between 3 points might include the equivalent of U-Turn in the middle of the street which on in reality is not possible). Did anybody solve this kind of problem, preferably with opensource tools (I know that ArcGIS have this kind of capabilities and I am trying to see if I can get there with it, but I generally prefer opensource)

Since these results are not optimum I end up using the graph road plugin (by the way is there a way to export a newly created path to an already existing file to create a full route, like it is possible in road graph with the network tools in QGIS 3?) but it is really time consuming.


texts: ArcGIS temporal relationship between two feature layers

texts:
how can I built a temporal relation between two feature layers.

My first layer has a observation time period and my second feature has the actual observation. I would like to build the relation directly in ArcGIS using the python console.

What I need is a function which:

iterates through the second feature and
checks if the time condition with the first feature attributes are fulfilled and
returns the corresponding id.

My features have the following structure:

First feature layer

ID    DATE        START    END
1     2015-03-29  13:20    14:50
2     2015-04-05  08:35    12:05

Second feature layer

ID    DATE        TIME    IDfirstLayer
1     2015-03-29  13:23   0
2     2015-03-29  13:28   0
3     2015-03-29  13:34   0
4     2015-03-29  13:49   0
5     2015-04-05  08:39   0
6     2015-04-05  10:04   0
7     2015-04-05  11:35   0

I guess, beside the loop function, I need to build a timestamp for each date.


texts: What is a valid way of calculating (bioclimatic) variables from climate data in terms of downscaling?

texts:
I have average monthly values of minimum and maximum temperature and precipitation and I want to calculate the bioclimatic variables analogue to WorldClim. I also do statistical downscaling with regression kriging using x,y and altitude from a digital elevation model. The question is:

Should I calculate the bioclimatic variables with the downscaled data or can I downscale the calculated bioclimatic variables afterwards?

What is the difference?


texts: Running out of memory with arcpy.da.UpdateCursor()?

texts:
Could anyone let me know if they see why this would cause a memory issue? The user running the script let me know he got a message that he ran out of memory while the script was running.

EDIT: Originally this was being used by someone else on the large dataset, but now I am running on my machine and getting an 'Out of Memory' message. No error message is given in the python window. I just get a small popup that says 'Out of Memory', or once it just shut down ArcMap. I watched the performace tab in the Task Manager and it the memory just slowly rises the whole time. It goes through ~5000 points before giving this message. I thought maybe deleting the variables after each row might help (just a guess), but didn't seem to.

This is being run with ArcMap in the Python window. It takes the layers from the table of contents as inputs. Not ideal, but needed to use SelectByDimension which updates the mosaic dataset in ArcMap.

Task: iterate through points, get a date, filter a mosaic dataset of rasters with this date, and then get the cell value at the point location. I am not familiar enough with these operations to know what is building up in the memory. Anyone have any ideas.

def getSST(mosaicds,detections,dateField,longField,latField,sstField,reportID):
count = 1
with arcpy.da.UpdateCursor(detections,[dateField,longField,latField,sstField,reportID]) as cursor:
    for row in cursor:
        split = row[0].split("/")
        timestamp = "{}-{}-{}".format(split[2],split[0],split[1])
        arcpy.SelectByDimension_md(mosaicds,[["StdTime",timestamp]],"BY_VALUE")
        coords = str(row[1])+" "+str(row[2])
        value = arcpy.GetCellValue_management(mosaicds,coords)
        SST = value.getOutput(0)
        if SST == "NoData":
            row[3] = -999
        else:
            row[3] = SST
        cursor.updateRow(row)
        if count % 500 == 0:
            print count
        count += 1
        del split, timestamp, coords, value, SST


texts: How to calculate area of polygon using ArcGIS system tools?

texts:
I like "Add Geometry Attributes" tool, it has some great options, e.g. centroid inside. Documentation states that:

This tool presents an easy and efficient way to automate the task that
can also be accomplished using the Calculate Geometry dialog box
or the Calculate Field geoprocessing tool.

Great news, unfortunately I cannot figure out how to calculate polygon area or length of a line.
When I start tool I can see them:

In documentation the description for selected couple says "AREA - POLY_AREA: The area of the polygon" and "LENGTH: The length of the line."
However, when I provide input as projected feature class containing polygons, these options disappear and I see only geodesic:

I know how to do it in script, field calculator but I'd like to use system tools because they are fast.
How do I calculate area of polygon using ArcGIS system tools?


texts: Tie line width to map scale by a formula?

texts:
I will explain it by an example :
In scale 8000 the width of a road is 2 mm = 16 map units.
If I use the unit mm, then width at any scale will be 2 mm.
In scale 32000 I find a width of 2 mm is too much.
If I instead use the unit map units, then the width at any scale will be 16 map unis.
In scale 32000 the width will be  16 map units = 0,5 mm, which I find too small.
In scale 32000 I want a width of 1 mm.
This I could get, if I could calculate the width by a formula :
width = sqrt(32000/scale) mm
Is that possible in QGIS ?


texts:
I have a project with many layers. Some layers are objects of interest while other layers are custom layers. Simplified example:

Streets (polygon)
Buildings (polygon)
Street names (point)

I created a print layout and now I want to add a dynamic legend. The legend should update automatically if I change the color of a layer, but I do not want to see some layers in the legends. I want only the polygons of the streets and houses, but not the street names point in the legend.

If I remove the check in "update automatically" in order to remove or add items the legend loses the update option. So if I change the color of the roads the legend shows the old color.

How can I update the legend with layer style changes, without losing my changes to which layers are displayed?


texts: Can I project 3D DXF vertical cross section into a "plan view projection plane" for display in ArcMap?

texts:
Below is the lenghty explanation but simply put my question is:
How can I project a 3D DXF file (vertical cross section) into a "plan view projection" for use in ArcMap?

I am using Target (by Geosoft) software to visualize drillhole data in order to produce cross-sections, 3D views, etc. of sub-surface exploration data. The software is capable of all kinds of things a traditional GIS such as ArcGIS simply cannot do but unfortunately the final output is not very polished. There is very limited control over what I can do with the layout, labels, graphics... Unfortunately, there isn't a way to export the finished cross section into a vector PDF which would allow me to tweak the output in Illustrator or CorelDRAW. (Using Print to PDF converts to image pdf)

I can only export the cross section into a 3D DXF which, when loaded into CorelDRAW for example, does not display the text annotations.

When I view the DXF in ArcScene it displays fine but it is "correctly" vertically projected and again lacking text annotations due to ArcScene limitations. When I open the file in ArcMap, the section does contain the annotations but it is "correctly" viewed in "plan view".

Bellow is the original section as it appears in Target as well as screen captures showing how the exported 3D DXF is handled byt the various other applications I have at my disposal (ArcScene, ArcMap and CorelDRAW) I also have Leapfrog Mining but have just installed it so I am not up to speed with this software yet but if anyone knows of a way to utilize Leapfrog please feel free to provide suggestions.

Original Section in Target:

Exported DXF in ArcScene:

Exported DXF in ArcMap (You are basically seeing the section plan view):

And finally DXF in CorelDraw (would be usable but not perfect as it would require a lot of manual work. and.. if it included text):


texts: creating a new vector layer from a part of the current one

texts: I am using qgis 2.6 for the first time. As I was using 2.2 previously, it was easy for me to create a new vector layer from a part of the current layer with the help of option "save selection as a new vector layer" from Layer menu. But I don't seem to find this in qgis 2.6.


texts: Selecting lines intersecting lines on two points in QGIS 3.14

texts:
I have two shapes of lines: one A with exploded red lines, the other B with closed black lines.
In a virtual layer I want to select the yellow ones (A) i.e. those having only 2 points in common with B.
I try to find a way to count the number of intersections between 2 lines without success.
Or maybe there is a better way to do it.

@JGH There is something weird when I work with my data. I can't get the result you show on your screenshot.    When I run the query to get the same result as you:
SELECT poly.*
FROM input1 as poly, input2 AS ouverture
WHERE st_intersects(poly.geometry, ouverture.geometry)
AND st_dimension(st_intersection(poly.geometry, ouverture.geometry)) = 0

Then the result is the red lines :

The green polygons (poly) have only 4 nodes and topology is correct with black lines.    It runs from a model so green polygon file converted to lines is not visible.
The strangest are the two red lines at the bottom which are the opposite of the expected result.
@Taras : result is empty.
EDIT : sample of data
shp and comments inside.


texts: Far side general perspective / custom orthographic projection in QGIS?

texts:
I'd like to create a map with a far side general perspective as shown on the right:

Unsure if relevant but I'd also like to orient the perspective from GPS coordinates and zoom in, creating a realistic view of a continent without compressing edges as much. Extremities of the projection / view aren't of concern.
I'm using QGIS 3.10 but can't find documentation on how to do this. Is this possible, and if so, how?


texts: OSM Carto Street Numbers style rules for OSM Bright

texts:
I'm running an OSM tile server stack using Mod_tile, renderd, mapnik, osm2pgsql and a postgresql/postgis database on Ubuntu Linux 14.04.

I've chosen OSM Bright for my tiles, but i would like to incorporate some of the style rules found in OSM Carto, just not all of it.

Specifically, i would like to incorporate the style for "Street Numbers", but i'm having trouble identifying which .mss stylesheet in OSM carto that this style can be found and copied.

In the attached images, i highlighted a street number which is rendered using the style rule in OSM Carto that i am trying to identify and copy into OSM Bright.

Which style rule is used to render these street numbers?


texts: glcm R package correlation metric values outside bounds of -1 to +1

texts:
I've been using the glcm package in R to calculate various texture measures of raster data sets; however, I've been getting odd values for the correlation measure which should be constrained to -1 to +1. Although most values fall in this range, some pixels are assigned -Inf or Inf, and some pixels are assigned real numbers greater than 1. I've searched but haven't found anything about this issue and was wondering if anyone knew what would cause this and how to account for it. Below is a small reproducible example using R's volcano data set:

library(raster)
library(glcm)

rast<- raster(volcano)
textures<- glcm(rast)

hist(textures$glcm_correlation) #Most values in valid range (-1 to +1)
cellStats(textures$glcm_correlation, max) #Maximum Value is Inf
#[1] Inf
cellStats(textures$glcm_correlation, min) #Minimum Value is -Inf
#[1] -Inf

head(sort(unique(textures$glcm_correlation@data@values), decreasing = TRUE)) #There are real number values >1
#[1]      Inf 5.839971 1.000000 1.000000 1.000000 1.000000
head(sort(unique(textures$glcm_correlation@data@values), decreasing = FALSE)) #-Inf is <0 but there are no real number values <0
#[1]       -Inf -0.5735393 -0.5000000 -0.5000000 -0.4588315 -0.4335550


texts:
I need a map that can represent density estimates. The estimates are tied to points, but are for 1 hectare (100m*100m) areas of which the point is the centre. Each point has at least one, but usually 3 or 4 other points that fall within its 1 ha area. I would like all of the overlapping areas to show up as means of the relevant estimates.

In other words I would like the density estimate tied to the point to be applied across 1 ha, and for the estimate to be averaged anywhere that it intersects the area of another estimate.

Is this possible?

Preferably I would be able to do this in QGIS, but am also open to instruction in ArcGIS.

So far, I have tried a point density map, and a kernal density map in arcgis, but the values were not averages of the estimated densities tied to each point.

I also tried converting the points to square polygons using the buffer tool and then the feature envelope to polygon tool in ArcGIS, which allowed me to make the points into 1 ha boxes. I was unsure what to do after that.

My question is also somewhat similar to this one from two years ago, which never received an answer: http://gis.stackexchange.com/questions/87429/how-to-simulate-venn-diagram-of-buffers-using-qgis


texts:
The Sinca WMTS service works great on QGIS 2.18 but you can not connect to it on QGIS 3.1.

http://gis.sinica.edu.tw/worldmap/wmts

When I try to connect to it on QGIS 3.2.2 I get an error saying "failed to download capabilities". When I connect to it on QGIS 2.18 it works perfectly

Is anyone else experiencing this and has anyone found a solution?

I have attached a screenshot of my network settings.


texts: Considering holes/constraints in Voronoi polygon creation in QGIS?

texts:
I am trying to create voronoi polygons in QGIS that would consider "holes" in the general domain. An example would be:

I actually created the Voronois in this image using QGIS through the GRASS command, then using the "Difference" tool to create the holes. A separate polygon shapefile, which contains the extents of the holes, was used as the "Difference" layer. An example application would be creating polygons around sampling points that were collected between structures that should be excluded from the analysis.

Two problems arise here:

The "difference" function does not appear to work 100% properly, with some polygon boundaries extending into the "holes". This can be fixed by finding row in the Attribute Table which do not have a polygon ID number (or ID of "0").
This type of after-the-fact "hole-punching" can result in discontinuous polygons, as shown by the red arrow in the image.

My question is: is there a Voronoi tool or plugin that can consider the presence of "holes" in the center of the domain, as a one-step process, and also eliminate the generation of discontinuous polygons? I envision that such a tool would extend a polygon boundary to the nearest intersection with another boundary, unless that other boundary slams against a "hole" boundary first.


texts: Referring to Null value with ArcPy Expressions?

texts:
How do I refer to Null values in a field using an expression in Python?

I've tried using the following as well as = None and =Null but nothing has worked yet:

ScoreCondition = "Cond_Score(!"+Condition+"!)"
codeblock1 = """def Cond_Score(cond):
                if cond > 90:
                    return 1
                elif cond >= 71 and cond <= 90:
                    return 2
                elif cond >= 51 and cond <= 70:
                    return 3
                elif cond >= 20 and cond <= 50:
                    return 4
                elif cond < 20:
                    return 5
                elif cond is None:
                    return 1"""


texts:
I want to delineate a watershed using ArcSWAT tool. Upon uploading the DEM from disk a window pops up & says:
Dataset must have a projected coordinate system. The current coordinate system is geographic. Please define a projected coordinate system of your DEM dataset using ArcToolbox before proceeding.
Where as the coordinate system of the DEM is GCS_WGS_1984. Could you help me i fixing it?


texts: Get id of last added feature in qgis python action

texts:
I created a python qgis action that open the dialog form so the user can add a new feature.

maintenance_lyr = QgsMapLayerRegistry.instance().mapLayer('od_maintenance_event20180212171049766')

maintenance_lyr.startEditing()
qgis.utils.iface.setActiveLayer(maintenance_lyr)
qgis.utils.iface.actionAddFeature().trigger()

My question: How can I get the id of the newly created feature?


texts: Points and polylines representation of coordinates in QGIS?

texts:
I am working on the "Jakarta Bay", Indonesia. I have a shapefile for the Indonesian coastline and I clipped my region of interest using an appropriate polygon. I also have points that are located within the bay.

When I want to represent both layers together, I realize it is not possible because they have different coordinates format (see images). I don't understand it as I chose manually to set the SCR to EPSG:4326 (OTF). The project SCR is the same.

If I make a zoom to layer I can see each layers individually, but it looks like they are not in the same location.

Looking to the coordinates show different results (see images).

I would like to have coordinates as for the points for both layers (East, North) because I can show them in the OpenLayers plugin.

I have the feeling my scale is a bit strange for the coastline layer...


texts: HTML5 maps with dynamic data instead of Flash

texts: I was wanting to create maps on our site similar to http://globaledge.msu.edu/Countries/Rankings except that we want the maps to be HTML/HTML5 instead of Flash. (The above site is the site that I'm working on also) The data is dynamic and gets updated periodically, but there's no set time, I just update the database. Any help would be much appreciated. I've looked at MapBox, but they seem to be one time only type of maps and I'd have to create a different map for each variable.


texts: Import of print compositions in QGIS 3.0

texts:
I installed QGIS 3.0 last week and would like to edit my projects with it. However when I try to import my saved print compositions the following message appears (see picture below).
Can you tell my why that is or how I can save my compositions in QGIS 2.18 so I can use them in QGIS 3.0?


texts: QGIS forms : field display based on the value of another field

texts:
I have a form designed in QGIS 3.8.1. The form is associated with a layer for inputting data. Many of the fields on the form are only relevant if the first field is equal to True.

Example:

Field 1: Is the location positive for mosquito larvae? True/False

Field 2: What stage are the larvae at? L1, L2, L3, Pupae

I do not want Field 2 to be shown unless the value in Field 1 is equal to True. How can I achieve that the Field 2 is displayed only when Field 1 is equal to True?


texts:
I have a series of shape files that I am trying to convert and enter into a postgis database so that I can pull the data out as GeoJSON but I am having troubles getting the data to enter with the correct coordinate values.

My source files say they are:

Datum:  NAD 83, CSRS98
Projection: UTM Extended Zone 13

And when I run gdalsrsinfo against the .prj file, I get the following:

PROJ.4 : '+proj=utm +zone=13 +ellps=GRS80 +units=m +no_defs '

OGC WKT :
PROJCS["NAD_1983_CSRS98_UTM_Zone_13N",
GEOGCS["GCS_North_American_1983_CSRS98",
    DATUM["NAD83_Canadian_Spatial_Reference_System",
        SPHEROID["GRS_1980",6378137.0,298.257222101]],
    PRIMEM["Greenwich",0.0],
    UNIT["Degree",0.0174532925199433]],
PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",500000.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",-105.0],
PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0.0],
    UNIT["Meter",1.0]]

Using that information, I pulled up http://georepository.com/crs_4140/NAD83-CSRS98.html and see the ID is listed at 4140.

So I ran

shp2pgsql -s 4140 shapefile > output.sql

But when I import the result into postgis, my coordinates are all out of sorts:

SELECT ST_AsGeoJSON(geom) FROM shapefile_import;
{"type":"MultiPolygon","coordinates":[[[[535472.058,5466281.22],[535069.83,5466278.55],[535068.17,5466514.07],[535068.03,5466534.18],[535066.992,5466680.332],[535064.14,5467082.11],[535061.3,5467483.903],[535058.46,5467885.69],[535460.67,5467887.69],[535862.87,5467889.69],[536265.07,5467891.69],[536667.28,5467893.69],[536669.945,5467492.576],[536672.61,5467091.46],[536675.278,5466690.345],[536677.54,5466350.32],[536677.81,5466310.46],[536677.95,5466289.23],[536276.118,5466286.56],[535874.28,5466283.89],[535472.058,5466281.22]]]]}

Can anyone point me in the right direction?  I have been using http://geojsonlint.com to validate.

Edit 1

JGH's answer is correct, I was using the wrong ID in my conversion.  It should have been 2151.

I ran now run SELECT ST_AsGeoJSON(ST_Transform(geom,4326)) FROM shapefile_import;


texts: Exporting map view of all added layers and basemap from app side using Google Earth Engine

texts:
Is there a way to have a download button link on an EE application that downloads not just an ee image in the viewport of an app, but all of layers as shown in the viewport and the basemap that is shown in the viewport?
The code described at Exporting image from app-side using Google Earth Engine  works, but exports only a single image object, but not other layers and not the basemap.


texts: Adding non georeferenced image as BaseLayer

texts:
I am trying to add non georeferenced image as a Base layer in OpenLayers3 (v.3.8.2). On top of that there is WMS, called from GeoServer, as Overlay.
I successfully loaded both layers, but they are not overlapping properly (there is a slight difference between them). Is it possible to do it without georeferencing static image? Also, do I need to state CRS for WMS in GeoServer when I am uploading it.

Here is my code (unfortunatelly, I am not able to give you a valid geoserver link):

 // projection EPSG: 3765
var projection = new ol.proj.Projection({
    code: 'EPSG:3765',
    extent: [250515.0793, 4698849.3024, 747014.5638, 5163391.4419]
});

// extent or bbox of layer in EPGS: 3765
var extent = [260515.0793, 4798849.3024, 757014.5638, 5263391.4419]

// baselayer: map
var imgWidth = 1025;
var imgHeight = 1024;
var url = 'path/to/the/image.png'

var source = new ol.source.ImageStatic({
    url: url,
    imageSize: [imgWidth, imgHeight],
    projection: projection,
    imageExtent: projection.getExtent()
});

var map = new ol.layer.Image({
    source: source
});

// overlay: layer
var layer = new ol.layer.Tile({
        layer:'layer',
    preload: Infinity,
    visible: true,
        source: new ol.source.TileWMS({
            url: 'http://localhost:8080/geoserver/workspace/wms',
            params: {'VERSION': '1.3.0',
                     'LAYERS': 'workspace:layer',
                     'TILED': true,
                     'bbox': 'extent',
                     //'WIDTH': 256,
                     //'HEIGHT': 256,
                     'CRS': 'EPSG:3765', //do I need to state here CRS?
             'FORMAT': 'image/png8'},
           serverType: 'geoserver',
          })
});

// view
var olview = new ol.View({
    projection: projection,
    center: [502024.4795, 4964030.8880],
    extent: extent,
    zoom: 2
})

// map
var map = new ol.Map({
    ol3Logo: false,
    target: 'map',
    renderer: 'canvas',
    view: olview,
    layers: [map, layer],
    controls: ol.control.defaults({}, [
    new ol.control.ScaleLine({
         units: 'metric',
    })
  ]),
});


texts:
I have my shapefile imported as shown below:

I calculate all all_simple_paths from root to all leaf nodes of a graph.
Red dot is a node to start with. It gives me all path from red dot to green dot.
But all_simple_paths do not return paths from node red dot to other leaves (in black). It comes blank.

Each line has start and end point geometry linestring.
Is there any condition in networkx to return all paths?

Code is as follows:

import networkx as nx
G = nx.read_shp(str(iface.mapCanvas().currentLayer().source()))
paths = []
for node in (x for x in G.nodes()):
  for leaf in (x for x in G.nodes_iter() if G.out_degree(x)==0):
    print node
    print leaf
    paths = list(nx.all_simple_paths(G, node, leaf))
    print paths


texts: Automatically setting standard Raster style when loading raster sets in QGIS?

texts:
For my project I work with a lot of raster data sets and when loading the data into QGIS they render with a white background which I would like to see transparent (no-data).
I also worked before with ArcMap and that had an option to set a standard Raster Layer rendering. Since I have hundreds of raster-files it is a bit tedious to set the style/transparency one-by-one.
Does anyone knows if QGIS has something similar?


texts: Counting attributes from line intersections with polygons on QGIS?

texts:
I have a line layer with 1300 features and a polygon layer with 19,000 features that represents buffers of points.

I need to determine how many polygons does each individual line intersect. And have a column in attributes with "buffer intersections" and another one with sum of attribute X

Like in the below example

The line in the right intersects 3 circles so it would have a value of 3 while the one on the left would have a higher value due to intersecting a lot of circles. It was possible to get this count by using the Combine tool in MMQGIS.

Let's say each circle has a value between 10-100 of attribute X, how to get the sum of this value in each line?

Is there any tool that does this or some other way by using QGIS?


texts:
I'm wondering if there's anyway to use a wildcard character in a field name. Like using ["Species_%" = 'Hemlock'] instead of ["Species_1" = 'Hemlock' OR "Species_2" = 'Hemlock' ...].

I tried the way that I described but it didn't work. Does anyone know if it's possible with another method within the Select By Attributes window?


texts: ArcGIS and .Net compatibility (Extending ArcObjects)

texts:
Backstory

Historically, when compiling a custom dll that extends ArcObjects you would compile on a PC with the lowest compatible version of ArcMap installed, for example 10.1.  Then that custom extension would work on ArcMap 10.1 and up until ESRI issues a new release that is a breaking change.  If you compiled on a higher version of ArcMap (say 10.5) and then tried to register it on a PC with a lower version of ArcMap (say 10.3.1) it would not work.

ArcMap 10.1 - 10.3.1 targets .Net 3.5 by default but can be changed to target .Net 4.0
ArcMap 10.4 and up targets .Net 4.5.

What puzzles me:

If I compile my custom dll extension with ArcMap 10.1 installed, my assembly has to target .Net 3.5.  But it still runs on ArcMap 10.4 and up which targets .Net 4.5.

Today, I compiled my custom extension targeting .Net 4.6.1 on ArcMap 10.5 and then installed it on an ArcMap 10.3.1 PC - expecting it to never even pass the ESRIRegasm.exe process - but it did and it all seems to work.

Question

Does anyone have a compatibility chart of some kind that shows what versions of ArcMap and .Net a custom extension that extends ArcObjects will work with including forward and backward compatibility of both?  I've looked and can only seem to find what each version on it's own supports, but that seems to contradict my findings above.

Below is a sample of what I'm looking for

Assembly targeting .Net 3.5 compiled on ArcMap 10.1 works on ArcMap 10.1 - 10.5.1
Assembly targeting .Net 4.0 compiled on ArcMap 10.1 works on ArcMap 10.1 - 10.5.1
Assembly targeting .Net 4.5 compiled on ArcMap 10.4 works on ArcMap ?
Assembly targeting .Net 4.6.1 compiled on ArcMap 10.5.1 works on ArcMap ?
and so on...


texts: Alternative statement for deprecated rendererV2().symbols()?

texts:
In my QGIS python plugin, this code triggers a deprecated warning (QgsFeatureRendererV2.symbols() is deprecated) :

myLayer.rendererV2().symbols()[0].changeSymbolLayer(0, mySymbolLayer)

What is the alternative up-to-date way to modify a layer's symbology ?
The doc says :

  use symbols( QgsRenderContext& context ) instead

but I can't find the correct way to follow that advice...


texts: How to quickly load two distinct areas into postGIS without using --append flag?

texts:
I have two distinct / non-overlapping areas that I wish to load into a postGIS database.

The trouble is that once one area is loaded, it becomes painfully slow to load the second area using the --append flag. (1.8k/s vs. the 200+k/s for the initial import using the --create flag.)

Is there a way to load two separate osm.pbf files into the same database tables without osm2psql checking for diffs in the process?


texts: QGIS 3x rotated label shifts direction without the reason

texts:
Could annyone help me with this problem.
I want to have number 22 in the same alignment as 23 is. Unfortunately when I hit 90 degree of rotation it shifts without the reason.
What is the reason behind it?


texts: Drawing the slope of linear regression in my chart

texts:
I can't converse fit to a list. It said the fit.reduceColums is not a function, and I can't even use the addBands to the fit or collection. I'm new to this and only want to add a slope line to the chart. I failed so many times to create the list. Link to complete code:
https://code.earthengine.google.com/5a78c37666a323bb90f5911551067901
An extract of functional code looks as follows:
var table = ee.FeatureCollection("users/tsukajoon/jxborderm84");

// link to the code that computes the Landsat LST
var LandsatLST = require('users/sofiaermida/landsat_smw_lst:modules/Landsat_LST.js')

// select region of interest, date range, and landsat satellite
// var geometry = ee.Geometry.Rectangle([-8.91, 40.0, -8.3, 40.4]);
var geometry = table.geometry();
var satellite = 'L7';
var date_start = '2001-01-01';
var date_end = '2021-01-01';
var use_ndvi = true;

// get landsat collection with added variables: NDVI, FVC, TPW, EM, LST
var LandsatColl = LandsatLST.collection(satellite, date_start, date_end, geometry, use_ndvi)
// print(LandsatColl)

var time1 = ['2001','2002','2003','2003','2004','2005','2006','2007','2008','2009'
            ,'2010','2011','2012','2013','2014','2015','2016','2017','2018','2019','2020'];
var time = ['2016','2017','2018','2019','2020'];
var list = ee.List([]);
Map.centerObject(geometry,6)
var cmap1 = ['blue', 'cyan', 'green', 'yellow', 'red'];
var cmap2 = ['F2F2F2','EFC2B3','ECB176','E9BD3A','E6E600','63C600','00A600'];

for(var index in time){
  var ymd = time[index]+'-01'+'-01'
  var image = LandsatColl.filterBounds(geometry)
                         .filterDate(ymd,ee.Date(ymd).advance(1,'year'))
                         .max()
                         .select('FVC')
                         .clip(geometry)
                         .setMulti({'system:index':ymd,'system:time_start':ee.Number.parse(time[index]),'FVCName':time[index]}
                         )
  // var image = image1.addBands(ee.Image.constant(1));
  // print(image)
  list = list.add(image)
  // print(list)
}

var imgYear = ee.ImageCollection.fromImages(list);
print(imgYear);

var createTimeBand = function(image) {
  return image.float().addBands(image.metadata('system:time_start'));}
// function createTimeBand(img) {
//   var year = img.date().difference(ee.Date('2001-01-01'), 'year');
//   return ee.Image(year).float().addBands(img);
// }
var collection = imgYear.map(createTimeBand);
var fit = collection.reduce(ee.Reducer.linearFit());
print(fit);
var fitinfo = fit.getInfo();
var scale = fitinfo.bands[0];
var s
print('scale',scale);
Map.addLayer(fit.select('scale'), {}, "FVCslope");

for (var i=0; i<5; i++) {
  var Time1 = ee.Image(list.get(i))
  // print(Time1)
  // Map.addLayer(Time1, {min:0.0, max:1.0, palette:cmap2}, "FVC"+i);
  Export.image.toDrive({
    image: Time1,
    folder:'JiangxiFVC',
    description: 'FVC'+i,
    scale: 500,
    region: geometry,
    crs:'EPSG:4326',
    maxPixels: 1e13,
    fileFormat: 'GeoTIFF'});
  }


texts: Graduated color display of overlapping polygons in ArcMap

texts:
Is there a way to display overlapping polygons in a graduated color symbology based on the level of overlap within the same feature layer? For example, in one area of the map there are two polygons that overlap, so the bottom polygon is light blue and the top polygon is medium blue, but in another area there is three polygons that overlap, the bottom is light blue, middle is medium blue and the top dark blue....etc in ArcMap
Basically the bottom most polygon is the lightest and the top most is the darkest.

EDIT-- Is this possible to do programmatically? ArcObjects?Python?


texts:
Possible Duplicate:
  How to get “.GDB” Esri File Geodatabase in Qgis for MAC OSX 10.6?

I have some data located here: https://nycopendata.socrata.com/Facilities-and-Structures/Building-Footprints/xe92-xce7

I'm trying to import this data into QGIS and am unable. I'm a novice at this. Ultimately, I'm trying to upload it to indiemapper.


texts:
I have a series of points around the greater london area (may be also around 100-200kms away from it) like those here in the format Lon, Lat:

43032.1,32681
95816.5,34473.8
15439.3,-13987

The Application to which those points are transferred uses leaflet but does not use any of EPSG3857, EPSG4326, EPSG3395 from http://leafletjs.com/reference.html unfortunately. At least that's what I found out when using the proj cmdline tool like this:

echo -48884.4 64535.7 | cs2cs -f "%.10f" +init=epsg:3395 +to +init=epsg:4326

I really tried a lot of possible projections and flipped the order but I cannot find a way to get a plausible result which should be something like 51.xxxx -0.1xxxxx

Does someone have any idea on what projection may be used here?


texts: Adding number of months (as integer) to date format in QGIS

texts:
I am creating a form for QFIELD in a QGIS project.
I have three fields:

Date of today ([$now]), as a date format,
Number of months allocated for a certain task, as an integer
The due date for this task, as a date format too.

What expression should I write into the default box to compute the due date? I have tried many things without any success.


texts: gdal and google earth (and my gps) not reporting the same value

texts:
I am playing with srtm files and I dont quite understand my results.  I am testing a lat lon of N39 W77 or 39.0,-77.0.  I have tested using the python script from this post: How to extract elevation from .HGT file? and using gdallocationinfo.  When I plot the point in google earth, it says the elevation should be 73-74m, gdal reports 81m.  I'm curious why the difference?

gdallocationinfo N39W077.hgt -wgs84 -77 39
Report:
  Location: (0P,1200L)
  Band 1:
    Value: 81


texts: Shoreline identification using ERDAS 2011

texts: How can we separate shoreline from shallow water sea and land area in a LISS 3 image using ERDAS 2011 software?


texts:
What is the difference between WGS84 and EPSG 4326?

It seems like for a given dataset it might be both WGS84 and EPSG:4326.


texts: Python API For Generating GeoJSON or WKT From Interactive Map

texts: I am looking for a Python API for generating GeoJSON or WKT interactively using a GUI based map interface. I have fun some website where one can generate the WKT polygon using by drawing a polygon over the area of interest but I would like to do this within my python window. Looking forward.


texts: Zooming rectangle not displayed in OpenLayers 3?

texts:
In ol 3.10, I am trying to use this code to zoom, but it is not displaying my zooming rectangle.

ol.interaction.DragZoom.dragZoom= new ol.interaction.DragZoom ({
  condition: ol.events.condition.always,
});

What is the possible reason?


texts: How to automatically fill a lat/long in attribute table

texts:
I need to plot the distribution of a plant species on a map. I have created a shapefile for the species plots (name of file is PLANT). The coordinate system of this file is WGS84. On the attribute table of this shapefile I put the fields Latitude, Longitude, ID and Type.

I have a georerefenced map of Madeira Island, using the WGS84 coordinate system. What I want to do is using the Madeira Island map as a reference, I want to plot the distribution of a plant species so that when I click a on the map to add a point to the PLANT shapefile, the lat/long coordinates of that point go automatically to the latitude and longitude fields on the PLANT attribute table.

Can anyone tell me how this can be done or reference me to a tutorial that explains how this can be done?

In short, how can I make the latitude and longitude fields of a attribute table be automatically filled when I click with the mouse on a map?

I'm not experienced in using QGIS.


texts: MMQGIS Geocoding gives Failure connecting to API: HTTP Error 400: Bad Request

texts:
I used Google with API to geocode my address. Everything is fine until yesterday I keep having the following error message:
Failure connecting to API: HTTP Error 400: Bad Request
I double-check my CSV file and address format, everything is consistent.
What can I do to solve this problem?


texts: Having trouble concatenating with Arcade across fields and having No Data when all columns are empty

texts:
I have data that needs to be in a pop-up that is contained across several fields.  If the fields have no data, I would like the pop-up to say "No Aviation Resources".  Here is what I have so far, it works except for the "No Aviation Resources" part....
var i = 0;
var features = [];

function addvalue(feat) {
    if (!IsEmpty(feat)) {
        features[i++] = feat;
    }
}
//Add your values
addvalue($feature.AirAttackSearch);
addvalue($feature.AirtankerSearch);
addvalue($feature.LeadPlaneSearch);
addvalue($feature.ReconSearch);
addvalue($feature.SmokejumperAircraftSearch);
addvalue($feature.RotorWIngSearch);
addvalue($feature.UASSearch);

IIf(IsEmpty(features), "No Aviation Resources", Concatenate(features, ";"))


texts: Enclose polygons that are not overlapped and remove gaps between them?

texts:
I have to create a big polygon with the outer ring from some selected polygons that are close to each other like in the picture.

Is it somehow possible two fill the small gap between the polygons as shown in the other picture

I'm using this query:

 SELECT ST_ExteriorRing((ST_DUMP(ST_Union(the_geom))).geom) FROM parcel

to get the outer line but get 4 rows result.

What do I've to do to get only one row of the outer line without the gaps?