texts: Debugging ogr2ogr "AddGeometryColumn failed" and "Terminating translation prematurely after failed translation of layer"

texts:
As noted in the title, I'm getting the following error:

ERROR 1: AddGeometryColumn failed for layer pretty_polys, layer creation has
failed.
ERROR 1: Terminating translation prematurely after failed
translation of layer mybeautifulshapefile (use -skipfailures to
skip errors)

I'm just adding a shapefile to PostGIS using ogr2ogr (PG connection data is all
fake placeholders, obviously):

ogr2ogr \
    -f PostgreSQL \
        PG:"host='000.00.000.00' port='5432' user='nullislandpatriot' password='nullisland4eva' dbname='data_i_like'" \
    "data/mybeautifulshapefile.shp" \
    -nln pretty_polys \
    -nlt POLYGON

I've tried to identify many possible sources for the error, but I can't seem to
figure out what's wrong, and the error isn't descriptive enough.

So my question is, how can I drill down into this error to find out why AddGeometryColumn is failing?

Here are some of the things I have already tried:

using different files. I receive this error with any file, including the
natural earth countries shapefile
checking the permissions in PostgreSQL. I have read and write access to
the public schema that is being used. Though I should note that I am not
the superuser, and we had to give explicit permission to allow me to access
spatial_ref_sys, and geometry_columns.
editing the -nlt option to make sure I'm using the correct geometry type.
checking different arrangements of -a_srs, -s_srs, -t_srs, and not
using them at all.
lots of other things before I tried the natural earth data, like
PRECISION=NO, different database schemas, text encoding, ...

My assumption is that the error has something to do with permissions or
settings in PostgreSQL or the syntax for the connection info, since I get the
same error for well-known shapefiles like Natural Earth. I've tried the
connection info withou the interior single quotes.

I also tried CPL_DEBUG=ON, which shows me a bit more of what is working:

Shape: DBF Codepage = UTF-8 for
/Users/bgolder/Downloads/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp
Shape: Treating as encoding 'UTF-8'.
OGR:
OGROpen(/Users/bgolder/Downloads/ne_10m_admin_0_countries/ne_10m_admin_0_countries.shp/0x7f9a79d00720)
succeeded as ESRI Shapefile.
PG: DBName="'data_i_like'"
PG: PostgreSQL version string : 'PostgreSQL 9.1.9 on x86_64-unknown-linux-gnu,
compiled by gcc (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3, 64-bit'
PG: PostGIS version string : '2.0 USE_GEOS=1 USE_PROJ=1 USE_STATS=1'
OGR_PG_NOTICE: NOTICE:  CREATE TABLE will create implicit sequence
"countries_ogc_fid_seq" for serial column "countries.ogc_fid"

OGR_PG_NOTICE: NOTICE:  CREATE TABLE / PRIMARY KEY will create implicit index
"countries_pkey" for table "countries"

ERROR 1: AddGeometryColumn failed for layer countries, layer creation has
failed.
ERROR 1: Terminating translation prematurely after failed
translation of layer ne_10m_admin_0_countries (use -skipfailures to skip
errors)

GDAL: In GDALDestroy - unloading GDAL shared library.

Specs:

PostGIS 2.0, with PostgreSQL 9.1.9 on Ubuntu 12.04
GDAL 1.11.0 (on OS X 10.9.3, installed via Homebrew)


texts:
I have a data source, that sometimes returns geometries which are not properly cleaned e.g.

the pink polygon is only 0.02 cm wide and is always on the edge of the whole geometry which consists of the brown polygon and the yellow one.

Is there a programmatic way I could prolong the yellow/brown border and with that line slice the pink one and merge the resulting two pink polygons to their neighbor?

So far I had to do this manually in QGIS using the vertex tool, but I would like to do this automatically using JTS or Geotools, as it is a bit time consuming.


texts: Check if path along triangular shapefile possible

texts:
I have a problem I came across but don't know how to solve it. I have a triangular shapefile of a river with one attribute value being the water depth. The idea is to find out if fish are able to swim from one end of the river to the other end. The condition that is required is, that the fish always has a water depth of at least 30cm along the path.
My solution was to first delete all the fields with water depth of less than 30cm. The next step I wanted to do then was to check if there is a path possible possible from point A to point B without a gap along the path. however I did not find out how to do this. Any idea on how to solve this? Or any better way on how to approach this problem?


texts:
I am trying to create a soil-type map of an area from point data using QGIS.  I tried interpolating the data, since soil is continuous.  I did this by allocating a number to each point type and then using the interpolation function.  I later created contours from the output file to delineate the areas.  The problem that I am having is that the results from the interpolation and contours do not line up with the point data.  The point data and contours are in close proximity to each other but do not overlay accurately.

I would like the output to be as if though I manually drew lines that connected and delineate the different types of soils/point data.  Now the points lie outside of the delineated area.  Is there another technique to do this or what am I doing wrong?

I used the Raster Interpolation Plugin in QGIS Wien.  The interpolation method was IDW, Number of columns =300, Number of rows =300, Cellsize X = 0.00004, Cellsize Y =0.00006


texts: postgis vector raster not being displayed error

texts:
postgis 2.3 postgres 9.5 64 bit windows 10.

I loaded a raster layer correctly into postgis

raster2pgsql.exe -s 3424 -I -C -M R:\Data\LUCZ_2017\ConservationPriorityAreas\data\rasters\abid -F -t 100x100 havdeel3 | psql -U postgres -d raster -h localhost -p 5432

and did a simple conversion to polygon which works correctly.

drop table if exists polys;
create table polys as
select st_polygon(rast)::geometry(multipolygon,3424) geom from havdeel3;

ALTER TABLE polys ADD COLUMN id SERIAL PRIMARY KEY;

This vector layer gets properly displayed in QGIS.

NOW I am trying to do the simplest conversion back to raster

drop table if exists raster_poly;
create table raster_poly as
select id,ST_AsRaster(geom, 50,50,'4BUI') rast from polys;

The table correctly displays but when I try and display it in QGIS it gives me this error

2018-07-10T15:03:35 WARNING  :
2018-07-10T15:04:47 CRITICAL    Layer is not valid : The layer dbname='scrap' host=localhost port=5432 user='postgres' password='--' sslmode=disable key='r_table_catalog' srid=3424 type=Polygon table="public"."raster_columns" (extent) sql= is not a valid layer and can not be added to the map. Reason:
2018-07-10T15:05:16 WARNING  :
2018-07-10T15:11:18 WARNING  :
2018-07-10T15:11:25 INFO     : Connection to scrap was successful.
2018-07-10T15:12:37 WARNING  :
2018-07-10T15:13:41 WARNING  :
2018-07-10T15:17:27 WARNING  :
2018-07-10T15:20:58 WARNING  :

I have seen this tutorial http://geospatialelucubrations.blogspot.com/2014/05/a-guide-to-rasterization-of-vector.html but it is not so helpful.

Seems to be a similar question to this List of vector to one raster in PostGIS which was never answered.

I also have tried the answer to this question with my problem and I still get the same error ST_AsRaster for polygon

Reworked query that doesn't work

select id,ST_AsRaster(st_transform(geom,3424), 100,100,ARRAY['8BUI'],
                  ARRAY[118]) rast from polys

or this doesn't work either

select id,ST_AsRaster(st_buffer(geom,0), 100,100,ARRAY['8BUI'],
                  ARRAY[118]) rast from polys

But how does this work?

SELECT ST_AsRaster(ST_Buffer(ST_GeomFromText('POLYGON((-30 40, -20 30, -25 20, -23 10, -30 40))', 4326), 50),100,100,ARRAY['8BUI'],ARRAY[118]);


texts: Changing NAD27 to NAD83 in R, Error in proj4string

texts:
I am having trouble changing latitude/longitude data (from spreadsheet dataset) from NAD27 to NAD83. I am using R. The spreadsheet has 2 million coordinates so I can't upload the data for you all to work with. This is the code I have so far:

lat=c(dataset$Latitude[which(dataset$Projection=="NAD27")])
lon=c(dataset$Longitude[which(dataset$Projection=="NAD27")])
d=data.frame(lat, lon)

coordinates(d)=c("lat", "lon")
proj4string(d)=CRS("+init=epsg:4267")

CRS.new=CRS("+init=epsg:4269")
Nad27to83_projection=spTransform(d, CRS.new)

However, after the proj4string line, I receive the following error:

  Error in proj4string<-(*tmp*, value = new("CRS", projargs = > "+init=epsg:4267 +proj=longlat +datum=NAD27 +no_defs +ellps=clrk66 +nadgrids=@conus,@alaska,@ntv2_0.gsb,@ntv1_can.dat")) :
    Geographical CRS given to non-conformant data: -9.8732102099e+10  3.8108480000e+06 -1.0300835057e+02  3.7622300000e+05

In addition, after the spTransform line, I receive this error:

  Error in spTransform(d, CRS.new) :
    No transformation possible from NA reference system

I am stuck. Please let me know if anyone can help!


texts: Invalid Number Error in Google Earth Engine?

texts:
I wrote the below code for sentinel5 data, but returns Invalid Number error in monthly loop. It used to work very well but now shows error. How can I fix it?

code link: https://code.earthengine.google.com/7df888dfff449b1dbcb99e64d7ffc255
// study area

Map.centerObject(table);

Map.addLayer(table);

// sentinel-5 carbon monoxide

var carbon = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_CO')
.filterBounds(table)
.filterDate('2019-01-01','2019-02-01')
.select('CO_column_number_density')
.mean()
.clip(table);

Map.addLayer(carbon)

// loop monthly

var start = ee.Date.fromYMD(2018,11,01);

var months = ee.List.sequence(0,6);

var startDates = months.map(function(d){
  return start.advance(d,'month');
});

// image collection

var carbon_monthly = function(m){

  var start = ee.Date(m);
  var end = ee.Date(m).advance(1,'month');
  var range = ee.DateRange(start,end);

  var carbon = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_CO')
  .filterBounds(table)
  .filterDate(range)
  .select('CO_column_number_density')
  .map(function(img){return img.clip(table)});

  return carbon.mean();
  };

var carbon = startDates.map(carbon_monthly);

var carbon_stack = ee.ImageCollection(carbon)
.toBands();

print(carbon_stack);

Map.addLayer(carbon_stack);

Export.image.toDrive({
  image: carbon_stack,
  description: 'carbon_data',
  scale: 1000,
  region: table,
  maxPixels: 1e9
});


texts: Retrieving lat-lon coordinate from projdef?


texts: EPSG 26716 On The Fly Issues in in QGIS NAD27-WGS84

texts:
QGIS Version 2.18.7 and ArcGIS 10.2

Working on a Project in Northern Nicaragua and have come across some issues in projections and transformations that are a bit above my head.

Nicaragua government uses EPSG:26716 (NAD27 Zone 16N)

Field Data Collected using Garmin Handheld GPS64
Datum set to NAD27 Central, Ellipsoid: Clarke 1866

QGIS Project set to EPSG:26716, OTF transformation enabled

Within QGIS project I import a Google or Bing satellite image using Quickmapservices. Confirm the project is still in EPSG:26716. I then import a .csv of points from a legal landsurveyor performed in EPSG:32616. I specify the x/y and the coordinate system. The points plot exactly as they should in relation to the image. Now if I import a .csv of the same points collected in NAD27CENTRAL from a hand-held Garmin (supposedly EPSG:26716) and specify the coordinate system as EPSG:26716 they plot a consistent 112m south of where they should be on the image.

The difference in y coordinate value between WGS84 and NAD27 is -196m

I've included an image illustrating this shift using a nearby point in a courtyard. The red star is in the correct position relative to the image. it was collected in WGS84.

The blue star was collected in NAD27Central (EPSG26716)

If I save the blue star as a new shapefile specifying the coordinate system as EPSG:26716 the resulting nad27 coordinates in EPSG:26716 are 542063mE 1456236mN.

In Summary

EPSG32616 coordinate in field 542063mE  1456327mN
EPSG26716 coordinate in field 542053mE  1456131mN
EPSG26716 coordinate from EPSG32616 in QGIS 542063mE 1456236mN

Using ArcGIS 10.2 I can start a new project in EPSG:26716, import the survey points in 32616, specify the coordinate system then transform them using NAD_1972_To_WGS_1984_2 and the points will line up with those collected using the Garminm NAD27 projection, however, this does not line up with basemap imagery form ArcGIS Online.

How do I tackle this discrepancy?


texts: Changing highlight property of layer or table field using ArcPy?

texts:
I need to highlight a field of a layer. With ArcMap it is accessed from Fields property page of a layer. From there one can reorder, hide and highlight the fields.

I couldn't find a property or method to change the highlight property of a field.

The FieldInfo Object only allows to set the visibility of a layer. Is there a way to change the highlight property?


texts:
I am searching for tutorials for the Image Segmentation modules of SAGA (version 2.0.8).
The official wiki from the software doesn't contain any info on this modules, same thing happens in the references (http://www.saga-gis.org/en/about/references.html).

Does anyone know where I can find some basic tutorials with examples on how to use those modules?


texts: I'm new about GeoServer. I want to download preview layer and data of the layer from GeoServer, how can I do that?


texts: QGIS Singleband pseudocolor Continuous symbology

texts:
If you want to configure the symbology of a raster in QGIS using render type 'Singleband pseudocolor', you can choose out of 3 modes to generate values in the table (continuous, equal interval and quantile). Although, I do not find information on the logic behind these modes. More precisely I'd like to know how the 'Continuous' mode classifies the values. In fact, my gut feeling expects this to be the same as 'Equal Interval' but obviously it isn't...
Could someone please point me to some documentation or explain what this 'continuous' mode does exactly?


texts: How to fix transportation network connectivity problem?

texts:
I am using network analysis in ArcGIS to conduct a health-related analysis of limited resource country. the quality of the road network layer is not high, with many connectivity problems between the road lines.

Probably there is no way to obtain a higher quality layer. What is the recommended way of testing and fixing the transportation network connectivity?

I use ArcGIS and R mainly but happy to learn other tools if needed.


texts:
I have two polygons in PostGIS that have incongruities. One polygon (the yellow one) is legal county boundaries. The other polygon (blue) is essentially neighborhoods. 98% of the lines align perfectly. The other 2% don't align because the boundaries have changes but not been updated on the neighborhoods polygon.

My question: Using either PostGIS (preferred) or QGIS, how can I "snap" the neighborhood boundaries to the county boundaries? I've tried:

ST_Multi(
      ST_Intersection(neighborhoods.geom, county.geom)
      ) END AS geom

EDIT: This query "worked" in that it prevented neighborhoods from bleeding into incorrect counties, but it just removes the area that causes the bleed creating blank spaces inbetween boundaries.


texts: How to merge big and many polygons (SHP) in QGIS/Spatialite?

texts:
Background:
I have a bunch of (huge) SHP-files that I would like to merge to one big SHP-file that I can use for geo-processing (clip, etc.). Now I was able to merge a big part of the SHP-files to one big ~ 400 MB SHP-file (Note: This big SHP-file is the result of merging many smaller SHP-files). But this is still incomplete, I would need to merge that big file with a couple of other SHP-files. Unfortunately this does not work in QGIS, it keeps crashing (minidump). I suppose it's because of the file format (not able to handle that amount of data). So I thought maybe it would be possible to add these SHP-files to a Spatialite DB in QGIS and then merge them there somehow. But I don't know exactly how to continue. Maybe there is an easier way to achieve the result

Question:
Is there a way to merge these (SHP-)files in QGIS-Spatialite or just QGIS so I can use them with other Layers (SHP-files) for exampling for clipping? I know that theoretically I could use the SHP-files that I want to clip with other SHP-layers one by one, but it would be nicer to have all the data in one layer/file/DB and only clip once. All the SHP-files contain polygons.
Edit: Basically the goal is to be able to geoprocess the data comfortably (clip, difference, etc.).

Below is some information from the attribute table from two huge SHP-files that I would like to merge (Merging SHP in Qgis does not work, it crashes):

SHP-file 1: Has one row, file size is 700 kB
SHP-file 2:Around 35,000 rows, file size is ~ 400 MB

Important Edit 14.8.18:

So, now I loaded all the many and partly big SHP-files (or the data from them) into a Spatialite data base (or a *.sqlite file). I have them there, I guess they are called "tables" now? I still would like to join/merge/dissolve (what's the right word?) all these different layers/tables into one layer/table that I can use to run geoprocessing commands like clip or difference with another layer in QGIS (in this case a SHP-file layer). I could do it one by one and use the result layer for each layer/table, but the goal is to have it done in one round. I thought that's what all this fancy software is made for? If it helps all the tables/files/layers in my Spatialite DB contain a "geom" -column that says "MULTIPOLYGON" in each row.

I try to do something with st-Union as described for example here, but I can't get it to work. Can anybody help or give a hint?

Edit 15.8.18:
I was able to merge two layers in Spatialite that contain the same amount of columns with this SQL code (found it here):

CREATE TABLE merged AS
SELECT * FROM layer1
UNION ALL
SELECT * FROM layer2;

But I still have other layers that I want to merge and they have a different number of columns which creates this error: "SELECTs to the left and right of UNION ALL do not have the same number of result columns".

Also I'm not sure how to add the other layers, I tried this code and then got the error message aobve:

CREATE TABLE merged AS
SELECT * FROM layer1
UNION ALL
SELECT * FROM layer2
UNION ALL
Select * FROM layer3;

The final layer (merged from all the individual layers) can continue all columns. Important is that it can be used with QGIS' geo-processing tools with other (SHP)-layers.


texts: Script PyQGIS without Interface in Ubuntu

texts: please can you help me? I tried many times to do a script of pyqgis but without interface qgis. I use QGIS 3.6.3 and Ubuntu 18.04, does Somebody know do this?


texts:
I'm creating a map with some layers, each one has a label, the problem is that, Mapserver is displaying the labels of bottom layers, so it is mixing all the label.

Is it possible to hide label that is under others layers?


texts: Setting time animation by number of frames

texts:
I'm trying to use the ArcGIS 10 animation system to produce a time series of points by year. There are only 65 years in the dataset, and I want to make it so there is only one frame per year.
I have tried setting it to play a certain number of frames and frame duration, or by total duration, but when it exports, it alters the frame rate and adds a number of duplicate frames when I export to AVI or as sequential images.

Is there a way to export only a single frame for each year? The interval and units are set to 1 year.


texts:
I'm trying to rotate an arrow based on an attribute but the arrow is off by 90 degrees.

Is there a method for doing a mathematical expression on the attribute within the sld?

I've gone over the docs and this thread but have not had any luck.  It so far everything I have found pertains to arrows on lines and not arrows alone

this is how Im currently rotating

 <se:PointSymbolizer>
            <se:Graphic>
              <se:Mark>
                <se:WellKnownName>arrow</se:WellKnownName>
                <se:Fill>
                  <se:SvgParameter name="fill">#236de5</se:SvgParameter>
                </se:Fill>
                <se:Stroke>
                  <se:SvgParameter name="stroke">#000000</se:SvgParameter>
                </se:Stroke>
              </se:Mark>
              <se:Size>30</se:Size>
             <se:Rotation>
                <ogc:PropertyName>heading</ogc:PropertyName>
            </se:Rotation>
            </se:Graphic>
            </se:PointSymbolizer>

After reviewing Using simple Math in SLD file on Geoserver? along with reading the geoserver documentation again it is still unclear on where in the function you choose your operator.  In my case I would like to use addition.  The above question is for multiplication. I would be multiplying by 90 not adding 90, the documentation doe not have comparative examples of mathematical functions being used.  I do not see in the documentation where <ogc:Mul> is described. I can assume that <ogc:Add> would work but where is this supported with documentation.


texts: How to create a catchment area using LSOA and point file data

texts: I have to create a catchment area (based on LSOA) for the 19 underground stations within Camden borough and it has to be in a polygon vector file. I know I have to somehow group the SOAs but don’t know which geoprocessing tool to use.


texts: Saving python Rtree spatial index to file?

texts:
I would like to save an Rtree spatial index to file. I tried using pickle.dump, but when I loaded the .p file using pickle.load, the bounds were way off. Here's an example:

import rtree, pickle
from rtree import index
iidx=index.Index()
iidx.add(0,(1,2,3,4))
f=open('rtree2.p','wb')
pickle.dump(iidx,f)
f.close()

When I try opening the file:

f=open('rtree2.p','rb')
uudx=pickle.load(f)
f.close()
uudx.get_bounds()

I get:

[1.7976931348623157e+308, 1.7976931348623157e+308, -1.7976931348623157e+308, -1.7976931348623157e+308]


texts: Postgis ST_Distance query returning inaccurate results

texts:
I am using the following query to retrieve records from a postgresql DB

@services = Service.select("*, ST_Distance(services.lon_lat, ST_GeomFromText('POINT (lat lon)', 4326)) as st_distance").order("st_distance").limit(10)

services.lon_lat is a geography type.

After implementing a "get directions" aspect to my app I have noticed a inaccuracy in the distances returned from the above query and the directions in google maps.

It's important to note I did expect a slight discrepancy taking into account the roads to the destination may not be a straight line but it seems the inaccuracy is the other way, for example the 1st returned item in this query says its 216km away where google maps says 181km, the 10th says 247km compared to google maps 203km...the further the distance the larger the inaccuracy...the 46th item returned if I choose to return that many says 424km compared to 298km.


texts: Calculate employment density using LEHD Census Data in ArcGIS

texts:
I have data from the Census Bureau's LEHD files that record number of jobs within blocks groups. I want to convert the number of data in the LEHD  into a feature that I can use to calculate block density of the density of jobs in the block groups.
 How do I do that? I have read something about using the Float to Raster tool but I am not sure how to use it.


texts: Is it possible to list SpatialGridDataFrame files to be used in a loop?

texts:
I have a lot (>50) of SpatialGridDataFrame files that I need to turn into raster files, stack with a set of raster files I already have and then extract data from random points from each layer of all of these individual stacks. I have written the code to do this for each file at a time and it works fine but I was trying to use a loop so I could run this for all of the SpatialGridDataFrame files in one go rather than one at a time.

I wasn't successful at listing the SpatialGridDataFrame files together to use as an input into the loop and get the following (example just with a few of my original files);

Layer1 <- list(SpatialGridDataFrame1A,SpatialGridDataFrame2A,SpatialGridDataFrame3A,SpatialGridDataFrame4A)))
Layer2 <- list(SpatialGridDataFrame1B,SpatialGridDataFrame2B,SpatialGridDataFrame3B,SpatialGridDataFrame4B)))

But this doesn't work correctly and when I try to call up the lists they look like this;

> Layer1
[[1]]
class       : SpatialGridDataFrame
dimensions  : 11, 44, 484, 2  (nrow, ncol, ncell, nlayers)
resolution  : 0.01788854, 0.01788854  (x, y)
Error in .printSpatial(x) : object 'nf' not found

Is there a way of listing these files together successfully so that I can then run the following loop?

for (i in 1:4)
 {
  L1 <- raster(Layer1[i])
  L2 <- raster(Layer2[i])
  depth <- raster(Depth_int)

rasterratio<- (L1/L2)

stack_all <- stack(L1,L2,depth,rasterratio)
assign( paste("layers",i,sep="_"), stack_all)
}

Or maybe use some other function like lapply to achieve this?


texts: Two identical geography objects produce different STIntersects results in SQL Server 2019

texts:
Consider I have two POLYON geography object, @g1 and @g2. One from WKT and the other one converted from binary:
DECLARE @g1 geography = 'POLYGON ((-45.70072144031528 70.79588950876575,-45.70072144031528 -32.671894242015554,202.6781848096847 -32.671894242015554,202.6781848096847 70.79588950876575,-45.70072144031528 70.79588950876575))'
DECLARE @g2 geography = CONVERT(sys.geography,0xE610000001040500000038D88CDAEFB251404CE17A3DB1D946C0ECE969A1005640C04CE17A3DB1D946C0ECE969A1005640C0AC47A1B0B355694038D88CDAEFB25140AC47A1B0B355694038D88CDAEFB251404CE17A3DB1D946C001000000020000000001000000FFFFFFFF0000000003)

SQL Server reports they are equal if I execute the following query:
SELECT @g1.STEquals(@g2) G1EqualsG2

| G1EqualsG2 |
| ---------- |
| 1          |

Now consider I have the following LINESTRING object:
DECLARE @line geography
SET @line=CONVERT(geography,0xE610000001040700000022625E9DF79F3D40E89D984C229E5A402D625EED3DA03D403C9E98C81E9E5A400C625E2C86A03D40F39D08B8199E5A4091625E9CBBA03D4055A00808069E5A4033625EFCE2A03D400E9E0828F59D5A4083625E7C26A13D408C9F08F0F09D5A405A625EDC7AA13D409F9E0834E99D5A4001000000010000000001000000FFFFFFFF0000000002)

If I try to create intersections for @line with @g1 and @g2, I get two identical intersections:
SELECT @g1.STIntersection(@line) G1Intersection,
       @g2.STIntersection(@line) G2Intersection,
       @g2.STIntersection(@line).STEquals(@g1.STIntersection(@line)) IntersectionsEqual

| G1Intersection     | G2Intersection      | IntersectionsEqual |
| ------------------ | ------------------- | ------------------ |
| 0xE610000001040... | 0xE610000001040...  | 1                  |

But, if I check whether @line intersects with @g1 and @g2 using STIntersects, I get different result:
SELECT @line.STIntersects(@g1) IntersectsG1,
       @line.STIntersects(@g2) IntersectsG2

| InsersectsG1 | InsersectsG2 |
| ------------ | ------------ |
| 1            | 0            |

I'm totally confused here. Why does STIntersects produce different result for a LINESTRING with two identical POLYGONs, while they can produce identical intersections?
Can anyone shed some light on this?

Edit
You may have noticed that several longitude values in the code are greater that 180°. It's because I'm using EFCore+NTS to perform query from client, and NTS performs counter-clockwise check on POLYGONs before sending the query.
But it seems that NTS performs this check on a planar coordinate system. It only works when longitude value does not jump across 180° meridian.
In the above case, if I perform the query with orignal longitude value (-157.3218151903153°), NTS will throw an error complaining that my polygon is not counter-clockwise.
To work around this issue, I have to add 360° to longitude values on the right hand side.
This is supposed to work in SQL Server, and as you can see from the example above, it actually works if I create the polygon from WKT. So this comes back to my original question again, why does STIntersects produce different result for a LINESTRING with two identical POLYGONs?


texts: Export an intersection of features as a raster Google Earth Engine

texts:
I am using Google Earth Engine to export some clipped images I am interested, but I need to export the region which I am using to clip the images as a geoTIFF, in this case, the feature called protectedAreas, I can't find any examples of doing this. Below is the code I use to export images, so instead I would like to export the feature.

// Define start and end ee.Dates.
var startDate = ee.Date('2000-01-01');
var endDate = ee.Date('2010-01-01');

// Define the regional bounds
var region = ee.Geometry.Polygon(
        [[[139.81147500927887, -22.584034738321385],
          [139.81147500927887, -39.52984428265981],
          [155.19233438427887, -39.52984428265981],
          [155.19233438427887, -22.584034738321385]]], null, false);

// Subset the to a Country .
var congo = ee.Feature(
  region
);

// Load a FeatureCollection from a table dataset: 'RESOLVE' ecoregions.
var ecoregions = ee.FeatureCollection('RESOLVE/ECOREGIONS/2017');

// Subset o the bounds of the ecoregion feature
// and other criteria. Clip to the intersection with congo.
var protectedAreas = ecoregions
  .filter(ee.Filter.and(
    ee.Filter.bounds(congo.geometry()),
    ee.Filter.eq('BIOME_NUM', 4)             // 4 = TBMF AUstralia, 1=TMBF Amazonia
  ))
  .map(function(feat){
    return congo.intersection(feat,ee.ErrorMargin(1));
  });

Map.addLayer(protectedAreas, {}, 'Eco region ');

var fire = ee.ImageCollection('FIRMS')
             .filterDate(startDate, endDate);

//print(fire)
// Clip and add a date band
var clipToRegion = function(img) {
  var clipped = img.clip(protectedAreas);
  return clipped;
};

var fire_clipped =fire.map(clipToRegion)

var scale = fire_clipped.first().projection().nominalScale().getInfo();

Map.addLayer(fire_clipped.first(), {}, 'Fire clipped');

Export.image.toDrive({
    image: fire_clipped.first(),
    description: 'Fires',
    scale: scale,
    region: region,
    fileFormat: 'GeoTIFF',
    formatOptions: {
      cloudOptimized: true
    }
  });


texts:
Possible Duplicate:
  Sort layers in QGIS table of contents?

I write a qgis python plugin ,and i try to connect keybord keys to reoder layers in qgis .
Are there a method in qgis allows me to retrieve the ordered list of layers and another that set the order?


texts: Using nested if statements in Python Parser of ArcGIS Field Calculator?

texts:
The solution in VB for the base process is at If/then VBScript/Python code equivalent to SQL IN ('x','y',z') expression in ArcGIS Field Calculator? works great but I was trying to nest 'if' statements and can't.

So pre-logic code of

if [grid_code] = 1 then
f2= 2
end if
if [grid_code] = 2 then
f2=3.5
end if

...etc. but this doesn't work.

What am I doing wrong?

I tried the second 'if' before the 'end if' and also tried 'elif' but these give script errors.

Also I am seeing a warning message about "The calculated value is invalid for the row with objectID=1" and when you click click 'yes' to continue it works fine.


texts:
We have ArcGIS 10 installed on a nice clean Linux box.  We have developed scripts to run on the Windows desktop development machines, and modified them to run on Linux.

Now, if we run the scripts on the windows boxes, they 'perform' within operational tolerances and perform correctly.  However, once we run them on the Linux envronment, they seem to work less well.  I am led to believe that ArcPy is synchronous and that calls to various objects have to finish computing, before the next call toprocessing can take place.  It seems not.  A call to aggregate points to a polygon, then to buffer the polygon, see's an empty buffer returned, when the polygon exists.  There is some latency between the Dev box and the ArcSDE instance, but surely if this is truly synchronous, it would only move on if the work had completed?

Does anyone have any experience of this?


texts: Referencing an external editable CSV file in a spatialite base

texts:
I would like to join a CSV file attributes to multiple spatialite layers, but I would like to be able to edit the csv externally (export it from a different project on Excel or Numbers every time I make a change).

I've done this with SHP layers and it worked fine, as long as I re-opened the project to refresh the join.

I'm getting started with Spatialite layers instead of SHP mostly because I want the ability to create dynamic views with simple queries.
I'm wondering if it's possible to keep such link with a csv that is stored outside of the database, or if I will have to ditch Excel/Numbers altogether and only work from within the DB. I would like to keep these tools for many other functions that are not spatial-related, and for which I'm a lot more proficient than SQL/DB (zero experience).


texts:
I created a point shapefile in QGIS 2.01. Then, I changed the style of the points so they display the contents of the Acceleration field in the attribute table. Works exactly as I require.

Now I would like to apply the same style to a second field in the table, but when I select that field in the Style dialog box, all the symbols are reset to their original display characteristics (colour, size). Is there a way to retain the current style definitions, and simply apply them to a second field? The value ranges are identical for both fields, so the same symbology could apply.

I tried to save the style using the Save Style button. Tried both options (*.qml and *.sld), but neither of them created a file. The *.sld option generated an error message that I should check the folder permissions, but that does not make sense because I've been using that folder to store the shapefiles.


texts: Difference between two features in same polygon in QGIS

texts:
I have a polygon and some features are below another features. If I would have two different polygons I would use the "Difference" tool so it would take the overlapping feauters to the "surface".
Does anybody knows a method (if there´s any) on how to do a kind of difference tool but with different features in the same polygon?


texts:
How can I select a recreation site location based on the criterion that it should AT LEAST be 500 m from a river (shapefile) and AT MOST 15000 m from the same river using Spatial Analyst Extension?

-Umair


texts: How to split a Polyline Geometry at a certain point in ArcObjects?

texts:
I have a polyline geometry and wish to split it into two different geometries, based on a given point on that line, in ArcObjects.

I have tried to find out how to do this, but to no avail. I found that there are three methods which are sort of related to what I want.

IFeatureEdit.Split Method This seems to do exactly what I want, but it Works on a Feature, and not on a Geometry. I Do not want to edit existing features, but copy the geometry from an existing feature, split it in two parts, and then process those parts.
IPolycurve.SplitAtPoint Method I though that this method would be useful, but it just adds a new vertex at the point, and does not split the polyline
ITopologicalOperator.Cut Method This would also work, but when I use it, giving a perpendicular line as the input paramter, I am getting an error which says 'An Internal error has occured in the geometry system.'

Is there some way of achieving what I want to do?


texts: Leaflet v1 plugin: layer control for ungrouped, grouped, and nested-grouped layers

texts:
I need a Leaflet plugin which extends L.Control.Layers. I have the following requirements:

works with Leaflet 1.x
handles ungrouped layers alongside grouped ones
handles nested layer groups
does not rely on actual Leaflet LayerGroups, but can be built up independently (in other words, layers can be grouped in this control in a different structure to how the layers themselves are grouped in the Leaflet map)

I've tried:

https://github.com/ismyrnow/leaflet-groupedlayercontrol (doesn't seem to handle ungrouped layers)
https://github.com/stefanocudini/leaflet-panel-layers (can't get it to work at all)
https://github.com/robbiet480/leaflet-categorized-layers (also can't get this to work)

Can anyone suggest if any plugin does what I need?

BOUNTY ADDED
To qualify for the bounty, please link to or upload a fully working example using Leaflet 1.x, including in the custom layer control:

ungrouped layers
grouped layers
nested groups


texts: Export raster with its min max range from GEE 2

texts:
Following up with my previous question Export raster with its min max range from GEE. I have tried the code below:
1.adapted from this answer Using computed min/max image values to set min/max visualization parameters in GEE
Use the .evaluate() function to convert server-side objects to client-side objects. Include the computed min and max values in an ee.Dictionary object and applied the .evaluate() function to it. A client-side dictionary object (dict) is made available within the scope of the anonymous function, where the min and max values can then be referenced and set as visualization parameters.

var minMax = ee.Dictionary({
  minVal: stats.getNumber('NDVI_min'),
  maxVal: stats.getNumber('NDVI_max')
});

print(minMax);

minMax.evaluate(function(dict) {
  var vizParams = {
    min: dict.minVal,
    max: dict.maxVal,
    palette: ['green']
  };

Export.image.toDrive({
  image: maskImage.visualize(vizParams),
  description: '21072016UNDVI',
  scale: 10,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});
});

It is not working.
Then I tried this code:
2.adapted from this answer Any way to have variable min and max in one visualization parameter?.

var minMax = ee.Dictionary({
  minVal: stats.getNumber('NDVI_min'),
  maxVal: stats.getNumber('NDVI_max')
});

print(minMax);

minMax.evaluate(function(val) {
  var vizParams = {
    min: val.minVal,
    max: val.maxVal,
    palette: ['green']
  };
// Export a cloud-optimized GeoTIFF.
Export.image.toDrive({
  image: maskImage.visualize(vizParams),
  description: '21072016UNDVI',
  scale: 10,
  region: table,
  fileFormat: 'GeoTIFF',
  formatOptions: {
    cloudOptimized: true
  }
});
});

It is still not working....
Does it have something to do with the Nodata or null data?


texts: Error when reading Layer from .gpkg with readOGR into R and writing it to PostgreSQL

texts:
About the SessionInfo() I work with:

R version 3.4.2
rgdal version:  rgdal_1.2-13
sf version: sf_0.5-4

The problem:

Recently I started to work with Geopackages .gpkg in R. The goal is to read a layer into R, do simple selection and write it to PostgreSQL.

I found two options to load the data from .gpkg into R:

Example, which works well:

The data can be downloaded from this GIS-Browser -> Select "Naturgefahrenkarten" and choose format .gpkg

natrualHazards1 <- readOGR(dsn= "pathto.gpkg", "WB_SYN_GK_F") # readOGR
natrualHazards2 <- st_read(dsn= "pathto.gpkg", layer = "WB_SYN_GK_F") # st_read

To write it into PostgreSQL I chose the pgInsert from the rpostgis Library

pgInsert(con, name=c("public",paste(c("tablename"),collapse="")),
     naturalHazards1,
     geom = "geom", df.mode = FALSE,
     overwrite = TRUE,
     df.geom = c("geom",paste(c("(MultiPolygon,2056)"),collapse="")), geog = FALSE)

Example, which causes error:

The data can be downloaded from this GIS-Browser -> Select "Amtliche Vermessung in Farbe" and choose format .gpkg

landuse1 <- readOGR(dsn= fgdb, "Bodenbedeckung_BoFlaeche_Area") # readOGR

  Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv = use_iconv,  :
    Multiple incompatible geometries: wkbPolygon: 1549; NA: 1268

landuse2 <- st_read(dsn= "pathto.gpkg", layer = "Bodenbedeckung_BoFlaeche_Area) # st_read

conn = dbConnect(PostgreSQL(), dbname = "dbname", user= "postgres", password = "password", host = "localhost", port = 5432)
st_write_db(conn, landuse2 , drop=TRUE)

The error is caused while writing to the PostgreSQL db.

  Error in postgresqlExecStatement(conn, statement, ...) :
    RS-DBI driver: (could not Retrieve the result : FEHLER:  Spalte Â»geometrieÂ« existiert nicht
  LINE 1: ...bedeckung set wkb_geometry = ST_GeomFromEWKB(cast(Geometrie ...
                                                               ^
  HINT:  Vielleicht wurde beabsichtigt, auf die Spalte Â»bodenbedeckung.GeometrieÂ« zu verweisen.
  )

Hint:

In order to test the object of class sf "landuse" I tried to make a simple conversion to a SpatialDataFrame using this code:

landuse_sp <- as(landuse2, "Spatial")

  Error in as_Spatial(geom, IDs = row.names(from)) :
    conversion from feature type sfc_GEOMETRY to sp is not supported

Questions:

What's the problem of readOGR in the second example, how can I learn more about the problem in order to resolve it.
What do I have to consider writing objects of class sf to PostgreSQL?
Could the errors be caused from a storing problem?
Somebody has an idea how to read the second example into R and save it from there to PostgreSQL?


texts: Multiple Layers Boundless Suite-SDK ol3view

texts:
I am currently building a simple web application with the Boundless OpenGeoSuite SDK ol3view (Openlayer3 viewer) and I was able to get one layer to show up from my geoserver instance but I would like to know how to show more than just one layer.  This is the code in the app.js that I have so far of the layer showing right now. The other layers are in the same workspace on geoserver with varying geometry field types.

// ========= config section ================================================
var url = '/geoserver/ows?';
var featurePrefix = 'Campus'; //workspace
var featureType = 'Basemap'; //layer name
var featureNS = 'Campus1'; //namespace
var srsName = 'EPSG:900913'; //projection
var geometryName = 'the_geom'; //geometry type
var geometryType = 'MultiLineString'; //field type
var fields = ['layer']; //columns
var layerTitle = 'layer';
var infoFormat = 'application/vnd.ogc.gml/3.1.1'; // can also be 'text/html'
var center = [-10764594.758211, 4523072.3184791];
var zoom = 3;

// =========================================================================

Does anyone have an example in the code for the app.js that would show how to add multiple layers and or a group layer?


texts: Serving 4 band GeoTIFF in GeoServer 2.1.4

texts:
I have an uncompressed GeoTIFF file (figure 1) with the following band information:

Image Structure Metadata:
   INTERLEAVE=PIXEL
Band 1 Block=5000x1 Type=Byte, ColorInterp=Red
Band 2 Block=5000x1 Type=Byte, ColorInterp=Green
Band 3 Block=5000x1 Type=Byte, ColorInterp=Blue
Band 4 Block=5000x1 Type=Byte, ColorInterp=Undefined

Figure 1

I have compressed the file using the following set of parameters:

gdal_translate -of gtiff -co compress=jpeg -co jpeg_quality=85 -co tiled=yes -co photometric=rgb uncompressed.tif compressed.tif

This gives me a file with the following information:

Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined

Problem:
When I serve my compressed image through GeoServer 2.1.4 I get image in figure 2 instead of the image in figure 1, which is what I would like. My goal is to have GeoServer just display the first three bands as RGB and ignore the 4th band.

Figure 2

I have tried to style the raster layer in GeoServer using SLD by adding the following channel configuration, but it didn't help.

<ChannelSelection>
    <RedChannel>
        <SourceChannelName>1</SourceChannelName>
    </RedChannel>
    <GreenChannel>
        <SourceChannelName>2</SourceChannelName>
    </GreenChannel>
    <BlueChannel>
        <SourceChannelName>3</SourceChannelName>
    </BlueChannel>
</ChannelSelection>

Why is my image server through GeoServer 2.1.4 looking the way it is and how do I configure GeoServer to display the image 'correctly'?


texts:
I'm trying to calculate viewshed from a set a points (gpkg). From what I can understand, I first need to run the "create Viewpoints" algorithm. After I do it, a new layer named "output" is created, but I don't understand what my next step should be. If I try to simply run the viewshed command (from the visibility analysis plugin), it does not work, indicating that no viewpoints exist in that area.

What am I missing here? This used to be quite a simple task in versions 2.14 and 2.18.

As a workaround, I've tried manually creating observ_hgt and radius fields in the points layers and then runing the Viewshed command. It seemed to be working, but in the end I got the following error:

Traceback (most recent call last):
File "C:/Users/tiago/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\viewshed_raster.py", line 231, in processAlgorithm
interpolate = precision > 0)
File "C:/Users/tiago/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\modules\visibility.py", line 159, in viewshed_raster
dem.open_window (point["pix_coord"])
File "C:/Users/tiago/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\ViewshedAnalysis\algorithms\modules\Raster.py", line 347, in open_window
self.rst.ReadAsArray(*self.gdal_slice ).astype(float)
File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\site-packages\osgeo\gdal.py", line 1944, in ReadAsArray
callback_data = callback_data )
File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\site-packages\osgeo\gdal_array.py", line 263, in DatasetReadAsArray
callback_data = callback_data )
File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\site-packages\osgeo\gdal_array.py", line 339, in BandReadAsArray
buf_obj = numpy.empty([buf_ysize,buf_xsize], dtype = typecode)
ValueError: negative dimensions are not allowed


texts: New Feature Layer isn't showing up in Geodatabase

texts:
I've been using the following script to add a new layer to a geodatabase.

 arcpy.env.overwriteOutput = True
 in_features = "lease_polygon_testing"
 out_layer= "lease_polygon_testing_Layer3"
 where_clause="Created_On > CURRENT_DATE - 830"
 arcpy.MakeFeatureLayer_management(in_features,
                              out_layer,
                              where_clause,)
 arcpy.CopyFeatures_management(out_layer,
 "X:/UserFolders/AnneGilligan/WC_74_Testing/landtrac_lease_test.gdb")

The script returns the exit code 0, however, I cannot find the layer in the geodatabase.


texts:
I have the latest Qgis2Web and am running Qgis 2.2.0 - I was able to run the plugin and the programme went through the motions resulting in "success" (message box) and the resulting "folder".

I saved it to my local folder and opened the "index" file on both my chrome/safari - both resulted in a blank screen.

Is there something I have done wrong or missed out?

Here are some screen shots: (enclosed)


texts: How to select by attributes and calculate aggregates using Modeler?

texts:
I need to know how I can select/aggregate equal attributes in a table.

This is quite the same I want to do. But I want to do it with the modeler.  How to group and count attribute data?

The column "Number" has attributes like 1.1 or 1.3 or 4.3. There are a lot of 1.1 or 1.3....in the column. In step1 I want to find all the 1.1 or 1.3 ...for calculating later the weighted average of a value by using the area. The problem is that the "Numbers" are different in all my tables. In EXCEL I could use a data filter to select all the equal attributes. Its not easy to describe...sorry. To say it in one sentence: If „Numbers“ are equal than calculate the weighted average of the Values and create a new table

But hopefully the screenshot describes the problem better than my words.


texts: Selecting features on GeoNB layer in QGIS?

texts:
I have installed QGIS 3, I have uploaded layers from Service New brunswick (GEONB), they all work fine, I am specifically trying to select a land parcel within the "parcels" layer. I have selected the layer in the list and the "select features" buttons and "start editing" buttons are greyed out.

Any ideas how to make the editing or selecting tools work?


texts:
I have used xarray on some satellite data that I'd like to mask it using a shapefile. I have previously been doing this using a combination of PIL, numpy and gdal, but it would be cleaner to do this directly on the xarray dataset with Rasterio and GeoPandas.
I have the following xarray Dataset:

Which looks as follows:

I've opened my shapefile in GeoPandas using:
sf = geopandas.read_file('GRI_jergetal.shp')

I then want to mask my xarray dataset using the shapefile in GeoPandas, and it seems like rasterio.mask.mask has this capability. However, the code that I have used gives the following error (adding in a nodata and a transform attribute to my xarray as it did not contain those parameters originally).
sf = sf.to_crs('EPSG:32643')
ndvi.attrs['nodata'] = np.nan
ndvi.attrs['transform'] = good_ds.affine
rasterio.mask.mask(ndvi.isel(time=0), sf.loc[0, 'geometry'])

[Output]:
`TypeError: 'Polygon' object is not iterable`.

I've tried a method posted in a previous post Python mask NetCDF data using shapefile (xarray & GeoPandas) & Mask Rasterio raster with GeoPandas shapefile, but the method doesn't seem right for my present case. Have I made a mistake in invoking the rasterio.mask.mask() function or should I try a different method for getting this to work?


texts: Change google maps center & water colour

texts:
I have a website I've desigend and I would like to show a map under a contact form.  Unfortunately the contact form goes over my map.  I created it in Google My Maps.  I tried changing the lat & long in the map by manually typing over the coordinated.  Doing so moved the map to the right position but when I get the embed code it centres it again.

Is there a way to retain these changes?


texts:
I've done the ESRI training: Creating Python Toolboxes Using ArcGIS 10.1 but is there anything more in depth? Since VBA is no longer supported in ArcMap 10 and up I'm trying to convert a cross section tool's code from VBA into python. The tool is here: http://resources.arcgis.com/gallery/file/geoprocessing/details?entryID=C83CC388-1422-2418-7F10-B4D3DF5F1EE6
Or if anyone knows of another way to make cross sections in arcmap that would be great.


texts:
I have data that contains energy demand information from postcodes(sometimes single and sometimes a range)the column POSTCODE_V and POSTCODE_T mean POSTCODE_FROM and  POSTCODE_TO respectively. here's a sample of the data

example of the postcode data is as follows

I'm trying to spatially join this data to the postcode shapefile without losing information.

I tried a join based on postcode from and then to merge the data by generating an additional column with sequence of numbers with repeating numbers if its a range of postcodes and dissolving the data based on this column. This is not working very well as some postcodes do not have energy data. Can someone help me figure out the best way to deal with this? In addition, the postcode file is huge. It has 450000 records.


texts: Converting self-intersecting polygons into distinct polygons

texts:
How to convert self-intersecting polygons into distinctive, well behaving polygons?

My polygon inputs consist of self-intersecting polygons, which fails the st_isvalid test, ie: the below statement returns some self-intersection rings:

ogrinfo -sql "select geom, st_isvalidreason(geom) from out" intersectingpolygon.shp -q

How can I sanitize the problematic polygons in intersectingpolygon.shp by breaking them into a few distinctive, well-behaving polygons, using GDAL?

PS:

There is a similar question on GIS, but I don't find it useful because the answers seem to deal with how to avoid the bad inputs in the first place, and how to convert MultiPolygon into a few Polygons. Nothing about how to convert self-intersecting polygons into well-behaving ones.


texts: How can I export a list of street names from a selection in ArcMap?

texts:
I am trying to get a list of street names and ranges for a specific area of my map.  My intention is to make a street index of sorts from this data.  I would like to export this data into excel if possible.  Does anyone know if this can be done?  If so, I would appreciate any help you can give me.

I'm using ArcGIS 9.3

Thank you.


texts: Creating new layer with points from selection in different shape on QGIS

texts:
I'm importing a shapefile with an attribute table with Lat/Lon for each point.

Is there a way to create a new layer with the coordinates from this shape as points?


texts: Converting CAD 2D section into PolylineM(Z) using ArcGIS Desktop?

texts:
How do I convert a 2D section drawn in CAD into 3D polyline in ArcGIS. By 3D polyline I mean PolylineM(Z).

I have 300 sections. Drafting it manually seems a bit laborious.

Any Ideas?


texts: Retrieve all points within a given radius using quadkeys

texts:
I'm new to GIS concepts and I'm learning about quadkeys and usage.

I am inserting few points with its quadkey into DB, and want to retrieve all points within a given radius using quadkeys.

Assume I have a point (46,46) whose quad key on a zoom level 23 is,
12122022101301320022121 (using this link)

Now I want to search all the points in DB that are at a radius of 100km to a reference point, say (47,47)

My idea is, use the given distance (100km) to find a zoom level, form a quadkey for (47,47) with that level, and query the DB with this key to retrive all keys that has the given prefix.

Is this the right way of doing?
Or is there a better way to deal with this using quadkeys.


texts: What are the benefits of RGBI Geotiffs vs ECW format?

texts:
I just received a batch of imagery for my project area in 2 formats: Geotiff (RGBI) and ECW.  The imagery looks very similar and covers the same footprint, but the file size difference is massive: 270GB for the Geotiffs vs 11GB for the ECW.

Is there any added value in the RGBI Geotiffs I should be aware of?  Is it even worth taking up this kind of server space to store this data?

The intended uses include map making, landscape change monitoring, visualization, and supporting engineering work... but always looking for new & interesting ways to make use of the information we have.


texts: Adding ArcGIS MapServer Layer via Python

texts:
I'm trying to add an ArcGIS MapServer layer from python code. The code I am running (currently via the console to test) is:

data = QgsDataSourceURI()
data.setParam("url", "http://egis.fire.ca.gov/arcgis/rest/services/FRAP/DE_HHZTier2/MapServer")
data.setParam("layer", "0")

iface.addRasterLayer(str(data.encodedUri()), "mylayer", "arcgismapserver")

This results in the error:

Invalid Layer: AMSProvider Could not parse spatial reference Raster layer Provider is not valid (provider: arcgismapserver, URI: layer=0&url=http://egis.fire.ca.gov/arcgis/rest/services/FRAP/DE_HHZTier2/MapServer

Adding the same service/layer via the UI (Layer -> Add Layer -> Add ArcGIS MapServer Layer) works great.


texts: Custom arcade expression to change the colours of text labels on an ArcGIS Online web map?

texts: I have a web map where I am labelling the point features with the values in 2 different text fields (eg."id"-"operation type").  I'd like to change the colour of the text for the "id" field label to make it easier to distinguish on the map. Can this be done using an arcade expression?


texts: Detecting nearly overlapped polylines in ArcMap?

texts:
I have feature made of lines representing infrastructure of city (water, gas, electricity...).

They intersect, but in some areas they nearly overlap which causes some type of lines to be invisible in final print in paper in mesurement 1:5000.

How to detect nearly overlapped polylines in ArcMap?

I need a way to detect lines within e.g. 3m buffer distance from magenta line. Next step is to distribute detected lines within gray area that represent road so that lines can be visible in print.
In this picture below I would like to select part of dashed orange line because it's too close to magenta.


texts: Trace Network and use results for another new trace

texts:
Is it possible in arcpy to make different geometric traces?

Like when in trace A, a featurclass 1 is passed by the trace, a second trace B should begin at featureclass 1.


texts: Computing together "sum" and "if" formulas in GRASS Raster Calculator

texts:
I'm working on global climate projections raster maps. My starting data are global raster maps of daily future climate, where every map represent one single day (e.g. I have 365 maps for 1 year). Furthermore, for every cell in the single map I have a temperature value.
I want to create maps of indices of extreme climatic events (ECEs) based on my starting data. An example of ECE indices is "number of days over 25°C in a year" (SD, Summer Days).

I was able to do it for single maps. Use the "IF" operation in the Raster Calculator. The output is a binary global map for a single day, where 1 are the "t > 20°C" values of the map and 0 are the "t < 20°C" values.
I wasn't able to do it, not even with just two maps. "SUM" the binary maps all together (365 maps for 1 year, but possibly I would like to have a 20 year mean), so that I can get an indices as "number of days over 25°C in a year".

My real questions are:

Why the "SUM" formula does not work with just two raster (if the "IF" formula works)? (Less important question, I think I can solve it somehow)
Is it possible to compute together the "IF" and the "SUM" operations in the Raster Calculator? (Important question)
Is it possible to summarize together all the formulas related to one year (or more years) in just one calculation (for one year, it would mean: 365 "IF" plus 365 "SUM")? (Much Important question)


texts: OpenLayers2 - Introduction tutorial problem

texts:
I have just started looking into OpenLayers and I can't seem to even get the example provided to work. I just get a blank screen.

Am I missing something really obvious? Here's the HTML I am trying to run.

  <title>OpenLayers Example</title>


texts: List Esri GRID files from unique data sources in a CSV file using Python 2.7 and arcpy

texts:
I have a CSV file with a list of unique data sources extracted from an ArcGIS map document. I would like to list all ESRI GRID rasters that appear in the list.

One example is:

G:\GIS\Depts\Ex\PA\DB\Rasters\Grids\RDS\tzech_z

I understand the (arcgis.env.workspace) needs to be set first but I'm not sure how this can be done if there are multiple GRID files in different directories.

One code block I have tried results in a TypeError:

for pth in newList: # iterates through each unique path
  pthDir = os.path.dirname(pth)
  arcpy.env.workspace = pthDir
  esriGridList = arcpy.ListRasters("*", "GRID")
  print esriGridList

Produces the error:

  TypeError: 'NoneType' object is not iterable


texts: Why don't all columns appear when changing style of shapefile for cartography in QGIS?

texts:
QGIS 2.8.9

I have been working with this shapefile and have already successfully coloured it using a graduated style from my column 'TIV' - however recently when I want to change the attribute the shapefile is coloured by, only latitude and longitude are displayed as options when there are multiple columns in the attribute table, as shown in the images below.

I have tried reimporting the shapefile and still encounter the same problem.


texts:
While users using my application, they can view resource from browser and see my arcgis rest service URLs. So If I use token based security, security token can view.
Can I hide this from users?


texts: Calculate area, length, width of rectangles in shapefile with Python

texts:
I have a list of rectangles imported from a shapefile into Python in degrees. (EPSG:4326)

rectangles = [{'type': 'Polygon',
  'coordinates': [[(-61.55002968237319, 73.546895599812856),
    (-61.55006750140465, 73.54690126492163),
    (-61.55004657065852, 73.547002738155776),
    (-61.55000884141658, 73.546996653057328),
    (-61.55002968237319, 73.546895599812856)]]},
 {'type': 'Polygon',
  'coordinates': [[(-61.55006750140465, 73.54690126492163),
    (-61.55010523064657, 73.54690693003006),
    (-61.55008438977398, 73.547007961258084),
    (-61.55004657065852, 73.547002738155776),
    (-61.55006750140465, 73.54690126492163)]]},
 {'type': 'Polygon',
  'coordinates': [[(-61.55010523064657, 73.54690693003006),
    (-61.55014304972003, 73.546912595138153),
    (-61.55012211897393, 73.547013646160064),
    (-61.55008438977398, 73.547007961258084),
    (-61.55010523064657, 73.54690693003006)]]}]

What I would like to do is calculate the area of each rectangle in square meters, and then the length and width of each rectangle in meters.  There is no guarantee that the coordinates are in the same order for every rectangle.

What I would like to do is output a Python list for each metric, for example

area = [10.933836, 10.920366, 10.872225]
length = [4.721, 4.703, 4.725]
width = [2.316, 2.322, 2.301]

(Those are not the actual dimensions from the coordinates above, just an example of the output format.)

I tried using GeoPandas, but when I try to convert from EPSG:4326 to any other projection, the Python session hangs and doesn't complete, even for small files with only 3 shapes.  For example, I tried

import geopandas as gpd
test = gpd.read_file('test.shp')
print(test.crs)
> {'init': 'epsg:4326'}

test = test.to_crs({'proj':'cea'})

The last line never completes, even for small files, I'm not sure why.

What are some approaches I can use to calculate the area, length, width for multiple shapes and output to lists in Python?


texts: Persisting spatial data from pyspark df/rdd to postgis?


texts: Can arcpy.env.workspace be set multiple times in same script?

texts:
Is it possible to set arcpy.env.workspace multiple times in the same script to ListDataFeatures from different workspaces?

eg.

import arcpy

def readFromGdb():

arcpy.env.workspace = "C:/data/base.gdb"

gdbList = listFeatureClasses()

#do the work

def readFromShp():

arcpy.env.workspace = "C:/data_shp/"

shpList = listFeatureClasses()
#do the work


texts:
I am trying to find a way to calculate the direct problem on the WGS-84 ellipsoid but instead of using the initial Azimuth (e.g. lat1, lon1, distance, fwd azi), I want to compute it with the reverse Azimuth (e.g. lat1, lon1, distance, point2-reverse-azimuth).

Are there any tools/libraries out there to do that?

My purpose is to find the tangential latitudes and longitudes to a circle of radius X from any given point -- to calculate an accurate lat/long bounding box that contains all points within X distance.


texts: Getting a layer that is not active with PyQGIS

texts:
I want to retrieve attributes of selected objects on a layer with PyQGIS.
All example I see is for the active layer:

vlayer=qgis.utils.iface.activeLayer()
selected_bf = vlayer.selectedFeatures()

But, my features are not on the active layer. What is the way to say that vlayer is another layer but not the active layer?

I put this:

import processing
from qgis.core import *
vlayer = processing.getObject('mylayer')
selected_objects = vlayer.selectedFeatures()

as a filter expression in QGIS 2.8 but I can't open QGIS after.

Error: "AttributeError: 'NoneType' object has no attribute 'selectedFeatures'"


texts: from coordinates can i find from multilinestring from database

texts:
Hello i have the coordinates from two points(x,y) and i have a database with roads and i want to find with dijstrak algorithm wich road i take but the algorithm want source and target but in my application i take only (x,y) coordinates

SELECT seq, id1 AS node, id2 AS edge, cost
        FROM pgr_dijkstra(
                'SELECT id, source, target, cost FROM edge_table',
                7, 12, false, false
        );

Can sameone help me ?

UPDATE I try this but i take all gid from table roads and not the gid it is close of the point

SELECT gid, ST_Distance(ST_Transform(ST_GeomFromText('POINT(39.240443 26.208803)',3857),4326),geom) AS myLineDistance
FROM roads
WHERE geom && ST_expand(ST_Transform(ST_GeomFromText('POINT(39.240443 26.208803)',3857),4326), 100);


texts: When are experimental parts in Openlayers considered stable? What does it take: Do they need to be included in a number of versions, be included for a certain time period, be used by a given number of sites or perhaps be formally accepted by a part of the dev team?


texts: Getting all aliasnames of featureClassLayer in ArcObjects?

texts:
I'm a novice in programming:

I'm trying to create a code which gives me the alias name from the attribute table. This is the code below that I tried, but it gives me only the layernames instead of alias name:

Imports ESRI.ArcGIS.esriSystem
Imports ESRI.ArcGIS.Carto
Imports ESRI.ArcGIS.ArcMapUI
Imports ESRI.ArcGIS.Geodatabase
Imports ESRI.ArcGIS.DataSourcesGDB
Imports ESRI.ArcGIS.Display

Public Class Aliasname Inherits ESRI.ArcGIS.Desktop.AddIns.Button

    Public Sub New()
    End Sub

    Protected Overrides Sub OnClick()
        FeatureClassopen()
        AccessLayerViaMxDocument(My.ArcMap.Document)
        My.ArcMap.Application.CurrentTool = Nothing
    End Sub

    Public Sub FeatureClassopen()
        Dim wsf As IWorkspaceFactory = New FileGDBWorkspaceFactoryClass()
        Dim gdbFileAddress As String = "C:\worki\samples.gdb"
        Dim ws As IWorkspace = wsf.OpenFromFile(gdbFileAddress, My.ArcMap.Application.hWnd)

        Dim fws As IFeatureWorkspace = ws
        Dim fc As IFeatureClass = fws.OpenFeatureClass("abc")

        Dim fl As IFeatureLayer = New FeatureLayerClass()
        fl.Name = "abc"
        fl.FeatureClass = fc
        Dim mxdoc As IMxDocument = My.ArcMap.Application.Document
        mxdoc.AddLayer(fl)
        mxdoc.ActiveView.Refresh()
        mxdoc.UpdateContents()
    End Sub

    Shared Sub LayerViaMxDocument(ByVal pMxDocument As IMxDocument)
        Dim pMap As IMap
        Dim pMaps As IMaps = pMxDocument.Maps
        Dim i As Integer
        For i = 0 To pMaps.Count - 1
            pMap = pMaps.Item(i)
            Debug.WriteLine(pMap.Name)
            MsgBox(pMap.Name)
            Dim pEnumLayer = pMap.Layers(Nothing, True)
            pEnumLayer.Reset()
            Dim pLayer As ILayer = pEnumLayer.Next()
            Dim sb = New System.Text.StringBuilder()
            While pLayer IsNot Nothing
                Debug.WriteLine(pLayer.Name)
                sb.Append(pLayer.Name).Append(vbLf)
                pLayer = pEnumLayer.Next()
            End While
            MsgBox(sb.ToString())
        Next
    End Sub

    'Funktion not fuctional
    Function Numbreoftable(ByVal workspace As IFeatureWorkspace, ByVal tableName As String) As Integer
        ' table open
        Dim table As ITable = workspace.OpenTable(tableName)
        ' number of fields
        Return table.Fields.FieldCount
    End Function

    ' not functional
    Public Sub ShowDistinctFieldAliasNames(pFClass As ESRI.ArcGIS.Geodatabase.IFeatureClass)
        Try
            Dim fields As ESRI.ArcGIS.Geodatabase.IFields = pFClass.Fields
            Dim field As ESRI.ArcGIS.Geodatabase.IField

            For i As Integer = 0 To fields.FieldCount - 1
                field = fields.Field(i)
                If field.Name <> field.AliasName Then
                    Debug.WriteLine(field.Name & ", " & field.AliasName)
                End If
            Next
        Catch ex As Exception
            System.Windows.Forms.MessageBox.Show(ex.ToString, "Sub Error", Windows.Forms.MessageBoxButtons.OK, Windows.Forms.MessageBoxIcon.Error)
        End Try
    End Sub

    Protected Overrides Sub OnUpdate()
        Enabled = My.ArcMap.Application IsNot Nothing
    End Sub

End Class

I don't know what I'm doing wrong. Maybe someone can help me finding a solution.


texts: Calculating partial distances from whole segment in grid using QGIS?

texts:
I am using QGIS and trying to find a way to calculate partial distances of lines in a grid.

In the left part of the picture you can see the sum of the line lengths within the grid cell (108.96...), and the length of the whole yellow line/segment that intersect it in black (146.37...).

How do I trim the yellow line to multi-lines which will show the proportional length in each cell (i.e 108.96+...+Len(Xn) =146.37)?


texts: downloading data from UNdata onto ArcGIS 10.1

texts:
http://data.un.org/Search.aspx?q=population

I am trying to download data from UNdata and successfully opened one link of raster data, however the data that I really desire won't open on Arc. I hit the download button and it gives me the choices of XML, Comma, Semicolon, and Pipe. I can download the XML to a zipfile to the correct folder and drag it onto the Arc page but then it simply says "could not open specified file". It doesn't like the XML file. This data would be perfect for my project, is there any way I could get it to show up? or is it simply not compatible?


texts: GeoDjango/PostGIS distance lookup not working

texts:
I'm building an api to make distance lookups using Django 1.8, GeoDjango, Django REST Framework and PostGIS.

Locally I'm using PostGIS 2.1.7, my host (PythonAnywhere) uses 2.1.4... I'm not able to install an earlier version locally (not so strong with making and installing).

The problem is that locally an api request like http://localhost:8000/v1/test-entries/8.5513076999999988%2055.6684179999999955/ return nothing, as it should. But doing the same request (using the same coordinates) on the live version on PythonAnywhere return all entries in the DB. The code is the same both locally and on PA. The distance lookup code is like so:

return TestEntry.objects.filter(pnt__distance_lte=(pnt, 10))

The only difference I see is the PostGIS version, but it seems unlikely that such a basic functionality shouldn't work in an older version.

Does anyone have experience with these kinds of technologies and give me some hints to what could be wrong?


texts:
I have about/more than 400 GeoJSON layer files in my QGIS project. When I checked the task manager memory, it displays qgis.bin.exe *32 using about 1,334,000 of CPU memory. I had another version earlier which used shapefile instead of GeoJSON, and I remember the memory doesn't go beyond 1 million.

Will it be safer to use shapefiles instead of GeoJSON for all my layers?


texts:
Working with Leaflet Control Search. I would like to add a timeout function to the marker that is added to the map if a location is found.

Any ideas in which line I have to add "setTimeout" function?

I tried

showLocation: function(latlng, title) {

    if(this.options.zoom)
        this._map.setView(latlng, this.options.zoom);
        this._map.panTo(latlng);

    if(this._markerLoc)
    {
        this._markerLoc.setLatLng(latlng);
        this._markerLoc.setTitle(title);
        this._markerLoc.show();
        this._markerLoc.bindPopup(title);
        this._markerLoc = setTimeout(function() {
        this._markerLoc.hide()
        },1000);
        if(this.options.animateLocation)
            this._markerLoc.animate();

    }

    if(this.options.autoCollapse)
        this.collapse();
    return this;
}});

The result: marker still shows up...


texts: Using PyQGIS to create images from the QGIS print composer each showing a different layer?

texts:
I have many vector layers which I want to display and print, one at a time. The following script works well for producing an image from the map canvas, but I want to produce an image from a print composer and that part does not work - the images from the print composer show all the layers visible in the map canvas at the time the script is run.

from PyQt4.QtCore import *
from PyQt4.QtGui import *

outputFolder = "C:/Users/fsc/Documents/Temp/"

layers = iface.legendInterface().layers()
settings = iface.mapCanvas().mapSettings()

for lyr in layers:

  # display layer and wait for completion
  settings.setLayers([lyr.id()])
  job = QgsMapRendererParallelJob(settings)
  job.start()
  job.waitForFinished()

  # save image from map canvas
  image = job.renderedImage()
  image.save(outputFolder + lyr.name() + ".png")

  # generate image from print composer
  c = iface.activeComposers()[0].composition()
  c.refreshItems()
  dpi = c.printResolution()
  dpmm = dpi / 25.4
  width = int(dpmm * c.paperWidth())
  height = int(dpmm * c.paperHeight())

  # create output image and initialize it
  image = QImage(QSize(width, height), QImage.Format_ARGB32)
  image.setDotsPerMeterX(dpmm * 1000)
  image.setDotsPerMeterY(dpmm * 1000)
  image.fill(0)

  # render the composition
  imagePainter = QPainter(image)
  sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
  targetArea = QRectF(0, 0, width, height)
  c.render(imagePainter, targetArea, sourceArea)
  imagePainter.end()

  # save image from print composer
  image.save(outputFolder + lyr.name() + "-comp.png", "png")

I have tried turning the visibility of all layers off and then turning each on inside the loop before producing the output, but in that case the print composer just produces blank maps. (I think that would have worked before the introduction of multi-threaded rendering - MTR - but not now.) Any ideas anyone?


texts: Switching Costs in a Multi-Modal Network Database in ArcGIS

texts:
I am constructing a multi-modal transportation network using Network Analysis on ArcGIS. I have data on roads, railway, and rivers for the U.S. My goal is to calculate an Origin-Destination matrix of the cost to transport goods from each U.S. county to each port along the U.S.'s coastline.

What I would like to model is the fixed cost of switching modes of transportation. For example, goods transported using roads would incur a fixed cost of $500 to switch to railway transportation. This fixed cost may include the cost for unloading the truck, storage in a warehouse until departure, and loading them on a train.

What would be the simplest way to include this fixed switching cost in a work database?

I am using ArcGIS 10.2.2 for Desktop and the Network Analyst Extension.


texts: Questionable data on single channel satellite bands vs. RGB composites

texts:
I am using Sentinel-2 images for a given region. I have all the separate bands and a merged RGB file, which comes standard with a download from ESA's API.

Now I'm running some code where I separate data into two folders. This means tiny cutouts of the satellite image, which are distributed into either folder 1 or folder 2.

Now I can construct a histogram of f.x. Sentinel-2 Band 04(Red Channel), where two curves illustrate the mean values of the small images in folder 1 and 2. This is presented as such:

I've done this across all single channel bands (Except 02, which for some yet unknown reason wont work)

Then I've also run the entirety of the merged RGB file which looks as this:

So here is my question! Why does the single channel 04 band look so different, in comparison to the merged second picture, where the bands are separated in the code. FYI: Band 3 is also off.

Side info: The data I'm using is from 2016, which is before ESA included a standard merge. So the merged RGB file I'm using have been created through the Gdal script of merging exactly channel 4,3,2 together, which looks like this:

gdalbuildvrt -separate stack.vrt B04.jp2 BO3.jp2 BO2.jp2
gdal_translate stack.vrt merged.tif


texts: Can QGIS map canvas update while python script is running?

texts:
I want to automate a periodic refresh of the QGIS map canvas with updated shapefiles. A Python plugin or console script does not allow the map to update until the script terminates. Is there a way to force a map refresh in a running plugin or script, or is a stand-alone app required?

Using v1.7.4 now but could migrate to newer version.


texts: Is there maximum number of points you can import from CSV to QGIS 2.18

texts:
I tried to import a CSV file which contains 4,783 rows. However, only 20 records appeared. Any one has solutions to solve this situation? How can I add all rows in QGIS delimited layer?

PS: I use Mac 10.10.5, QGIS 2.18.0, GDAL/OGR 2.1.2.


texts:
I'm having an issue when offsetting a line shapefile. For some reason when the lines are offset, ArcGIS isn't drawing some segments although the segments are there confirmed via “identify”. Has someone come across a similar problem? why is this happening and how does it get resolved?

Thanks for the help


texts: Setting the 3d view function for a polygon layer

texts:
I am trying to create a 3d model of an archaeological site on QGIS. I had some contourlines, I made them polygons using the "polygonize" plugin and I categorized the polygons with a colour scale according to the altimetry value of each polygon.

At this point I used the 3d view in the layer properties panel to get a 3d view of the "polygons overlapping" effect. The value of the extrusion is the one found in the ELEV_2 field (i.e. the value of the height of each polygon).

The result is this:

There are two questions:

How do I keep the colour scale I used for the 2d symbology?
Is there a way to decrease the height of the base of the model?

I have already used Qgis2threejis. Unfortunately every time I try to open this layer with this plugin the pc freezes.
I have already tried this solution but applied to my case it does not respect the height increase of the polygons Jochen Schwarze's solution


texts: QGIS 2.12.3 Crashes on buffering line vectors

texts:
I'm using QGIS 2.12.3 and I'm trying to create a buffer of 100 mt around a line vector using the Vector - Geoprocessing Buffer tool, but qgis crashes at every try. I've already checked that the RS of the vector is the same of the project.
Can anybody give some hint to solve this problem and create the buffer?


texts: topographic correction using GRASS 7.0 (**i.topo.corr**)

texts:
I would like to apply topographic correction using GRASS 7.0 (i.topo.corr). The imagery were already converted from DNs to Top-of-Atmosphere Radiances. When I apply the i.topo.corr command it runs normally without any error, but the result is a "white" image. I cannot see the bands, it seems that I am missing something. I applied topographic correction for Landat images from 1991 and 1979 and it worked perfectly using the same commands, but for Landsat image from 2003 and 2015 it does not work.

I'll explain what I did step by step.

The illumination model was created as follows:
i.topo.corr -i --overwrite output=2003_111.84_SRTM_illum basemap=SRTM_1s_mosaic@PERMANENT zenith=54.56122270 azimuth=78.49871615 method=minnaert

After the illumination model is created, i used it to correct the desired bands of the image:
i.topo.corr --overwrite input=13Jan2003_111.84_toar_B1@PERMANENT,13Jan2003_111.84_toar_B2@PERMANENT,13Jan2003_111.84_toar_B3@PERMANENT, output=tcor basemap=2003_111.84_SRTM_illum@PERMANENT zenith=54.56122270 azimuth=78.49871615 method=minnaert


texts: python - psycopg2.errors.RaiseException find_srid() - could not find the corresponding SRID

texts:
I'm trying to import a layer into a PostGIS database with
gdf.to_postgis('test_polygons', engine, if_exists='append', index=False, dtype={'geom': Geometry(geometry_type='POLYGON', srid= 4326)})

And I'm getting the error:

Traceback (most recent call last):   File
"C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1276, in _execute_context
self.dialect.do_execute(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\default.py",
line 609, in do_execute
cursor.execute(statement, parameters) psycopg2.errors.RaiseException: find_srid() - could not find the
corresponding SRID - is the geometry registered in the
GEOMETRY_COLUMNS table?  Is there an uppercase/lowercase mismatch?
CONTEXT:  PL/pgSQL function find_srid(character varying,character
varying,character varying) line 17 at RAISE

I have a custom projection, which i set to the layer
Any ideas, what's wrong with geometry and srid?
Full error message:

Traceback (most recent call last):   File
"C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1276, in _execute_context
self.dialect.do_execute(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\default.py",
line 609, in do_execute
cursor.execute(statement, parameters) psycopg2.errors.RaiseException: find_srid() - could not find the
corresponding SRID - is the geometry registered in the
GEOMETRY_COLUMNS table?  Is there an uppercase/lowercase mismatch?
CONTEXT:  PL/pgSQL function find_srid(character varying,character
varying,character varying) line 17 at RAISE

The above exception was the direct cause of the following exception:

Traceback (most recent call last):   File "<pyshell#84>", line 1, in

gdf.to_postgis('test_polygons', engine, if_exists='append', index=False, dtype={'geom': Geometry(geometry_type='POLYGON', srid=
4326)})   File
"C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\geopandas\geodataframe.py",
line 1097, in to_postgis
geopandas.io.sql._write_postgis(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\geopandas\io\sql.py",
line 375, in _write_postgis
target_srid = connection.execute(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1003, in execute
return self.execute_text(object, multiparams, params)   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1172, in _execute_text
ret = self._execute_context(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1316, in _execute_context
self.handle_dbapi_exception(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1510, in handle_dbapi_exception
util.raise(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\util\compat.py",
line 182, in raise
raise exception   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\base.py",
line 1276, in _execute_context
self.dialect.do_execute(   File "C:\Users\Александр\AppData\Local\Programs\Python\Python38\lib\site-packages\sqlalchemy\engine\default.py",
line 609, in do_execute
cursor.execute(statement, parameters) sqlalchemy.exc.InternalError: (psycopg2.errors.RaiseException)
find_srid() - could not find the corresponding SRID - is the geometry
registered in the GEOMETRY_COLUMNS table?  Is there an
uppercase/lowercase mismatch? CONTEXT:  PL/pgSQL function
find_srid(character varying,character varying,character varying) line
17 at RAISE
[SQL: SELECT Find_SRID('public', 'test_polygons', 'geometry');]
(Background on this error at: http://sqlalche.me/e/13/2j85)


texts: How to reclassify raster data by their sign?

texts:
I want to do the following: I have two data raster sets with different units and values. One results from measured metric data, the other from modelled grid values. What I want to compare is simply whether the data sets differ in where they show positive and negative values.

So - all negative data on both sets should be "-1", all data between -0.001 and 0.001 should be "0", and all positive values simply "1".

Do you have an idea how to do this best on ArcGIS 10? I suspect raster calculator...?


texts: I have a GIS application where users have to login in with username and password. The WMS is served using GeoServer and displayed using OpenLayers. I need to restrict WMS layers.  I want to access the secured layers via WMS without exposing my credential. Is there any method for this? Just like URL filtering.


texts: Background Files disabled in Trimble GeoXT

texts: I am trying to send background shapefiles to the trimble GeoExplorer 2008 series. I have used GPS Pathfinder Office and sent the files in the correct coordinate system. Look at the picture I attached. The Geo7X clearly shows background files as an option. The GeoXT does not.


texts: Coverting text from one Column to a year date format in another column in QGIS

texts:
I have two columns called "LC" and "LC Datum".

"LC" contains text value's like 2015, 2016 and so on. These are actually years and need to be converted/transfered to another column called "LC DATUM", where the text value will be a date year value.

When I try to go in the Attribute table into the Field calculator. And then into the Column "LC Datum" and try to set the expression "LC" in the "LC Datum" Column it gives the value NULL. This is strange because I clearly see a result of the execution in the calculator window. I assume I need to do an extra conversion of my data. But how do I do this?

Edit: As the last thing I also need to note that there are null value's in my LC column.


texts: Carto, Named Maps, writing a list to check with 'in' sql format

texts:
I'm trying to use a Named Map in Carto with the 'in' operator to avoid to make various calls to paint the same way a few features.

I want to do this:

SELECT * FROM madrid where id_tramo in ('2805200000011','2805200000041','2805200000051','2805200000061')

In my case, I'm using a string with simple quotes and commas to separate values.This is part of the call:

  config: "provincia":"a_coruna","municipio":"15041","idtramos":"'1504100012461','1504100014191','1504100012951'","layer0":1}

You can see the parameter "idtramos" being "'1504100012461','1504100014191','1504100012951'". Of course, the Named Map is getting in trouble with simple quotes and commas because I'm using sql_literal as variable type. As we can see in Carto. And this is what I'm getting:

  errors:["syntax error at or near "1504100012461""]

It's possible to do it with a Named Map? If not I will have to make a loop calling the each feature, but I understand is not the best way to do it. Maybe there is another way to do it.

This is the Named Map:

'paint_bloqued_tramos' : """{
    "version": "0.0.1",
    "name": "paint_bloqued_tramos",
    "auth": {    ...    },
    "placeholders": {
        "provincia": {
            "type": "sql_literal",
            "default": "No hay variable provincia y no puede filtrar en la query del NM"
        },
        "municipio": {
            "type": "sql_literal",
            "default": "No hay variable municipio y por el NM da error"
        },
        "idtramos": {
            "type": "sql_literal",
            "default": "No hay variable idtramo y no puede filtrar en la query del NM"
        }
    },
    "layergroup": {
        "layers": [{
            "type": "cartodb",
            "options": {
                "cartocss_version": "2.1.1",
                "cartocss": "#layer { line-width: 7; line-color: #000000; line-opacity: 0.3; line-comp-op: color-burn; }",
                "sql": "SELECT * FROM <%= provincia %>_tram where municipio_id = '<%= municipio %>' AND id_tramo in (<%= idtramos %>)"
            }
        }]
    }
}


texts: How to compute the share of a line that is inside a polygon programatically?

texts:
I have two layers, one contains polygon while another consist of line. something like below figure:

Now I want to update line attribute column value through polygon attribute column only when, line is inside or contains or overlap 90% in the polygon. Like above picture it is depicted that polygon contains line almost 60%. How can I check that?

I only know that there are two function intersect and contains both cannot do this job. I tried it with contains but it didn't work as polygon provide area while line provide length so I can't compare both.


texts: Finding date and time of Digital Elevation Model

texts:
I downloaded a DEM from the Copernicus site.

I can provide the link: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1-0-and-derived-products/eu-dem-v1.0?tab=download

I would like to find the date of my DEM, so as to compute sun elevation and azimuth..

Do you have any suggestions on how to do that?

I tried to find a metadata file.
I was trying to find out an xml file or a text file, but when i found it, it doesn't include the specific information.


texts: Multiprocessing and ArcGIS in_memory workspace?

texts:
I am trying a script that makes use of multiprocessing and a number of tools within arcpy. The script works only if I write all input, intermediate, and output from child multiprocessing executions to the harddisk. Using the "in_memory" workspace as I had intended produces the feature class "does not exist or is not supported" error.

I've read over what I could find on problems with arcpy and multiprocessing. I've read several pages that specify that writing to the harddisk is necessary.

Is this true?

This seems like a horrible limitation of arcpy where combined with multiprocessing and I'm in disbelief that ESRI would let such a problem persist.

Can someone suggest a workaround?


texts: In which coordinates (x, y, z) are vertical grids expressed?

texts:
When working with 2D maps, coordinates are always expressed in X and Y. Adding a 3rd dimension also adds the Z-coordinate, which typically stands for height.

However, when working with vertical grids (for measurements, visualizations on building facades or the sides of 3D polygons, etc.) I don't know which coordinates to mention in order to express what I mean.

Is it 'X and Z' or 'Y and Z'? Or does it become 'X and Y' again since there is no 3rd dimension?


texts: ArcGIS Silverlight Web Mapping Application has two logins

texts:
I am trying to implement a Silverlight Web Mapping application with ArcGIS server and I'm trying to use the Token authentication with active directory as the user/roles store.

I have the service running so that I can successfully obtain tokens but when I try to initialize a map service the ArcGIS server challenges me again to authenticate with widows even though I have successfully created a token.

The IIS is set for windows authentication but in any other configuration (Anonymous, Impersonate) the map service never loads.

        //Import  Base Data Map Service
        ArcGISTiledMapServiceLayer BaseMap = new ArcGISTiledMapServiceLayer();
        BaseMap.Url = DataSourceStrings.ProductionServer + DataSourceStrings.I_GeneralWesternCanadaURL;
        BaseMap.Token = Token; // This is a valid token
        BaseMap.ID = "BaseMap";
        BaseMap.Visible = true;
        //BaseMap.ImageFormat = ArcGISDynamicMapServiceLayer.RestImageFormat.PNG32;
        BaseMap.Refresh();
        BaseMap.InitializationFailed += ArcGISDynamicMapServiceLayer_InitializationFailed;
        BaseMap.Initialized += ArcGISDynamicMapServiceLayer_Initialized;
        Map.Layers.Insert(0, BaseMap);

Solution...?
So I managed to fix this by setting the Server.dat file to disable security but the services still require a token in their web.config files. Seems Arcgis server like to handle the windows authentication but the design was to handle the authentication in silverlight using a wcf service. I've tested this and everything seems to work the way i want it too. Kinda of a hack way of doing it but everything works right.


texts: Using For Iterator in ArcGIS ModelBuilder?

texts:
I'm trying to build a rather easy model but somehow couldn't find a solution to one problem:
I work with a main model and a submodel which contains the iterator "for". I want to buffer a feature class ("building") with different distances. I added %value% to the output name so that the outputs are not overwritten. Afterwards I would like to use one of the outputs in the main model in a "union". But I have no idea how to get it as an input value. I thought I could open the tool and write down the path of the buffered feature class I want to use. But it's not working. I guess it's because it doesn't exists at that moment but I thought it'll work since it's created during the process of the model.

Is there any easy solution?


texts: vehicle tracking system in postgresql performance

texts:
I am developing the vehicle tracking system using postgresql.

When the data comes from device to my DB server, I am calling the following functions after inserting data into log table.

Distance Update -- picking up the last inserted row (latitude, longitude) of the specific device from log table and calculating distance with newly inserted row... updating the distance in the same table
Trip function -- gets the list of all locations'(stations) latitude, longitude, radius and again calculating the distance [by using this I am able to get the details of trip start date and trip end date and inserting those details into anotehr table]
Geofence functions
Overspeed function
some other functions

I will get almost 2000 rows per minute.

Will this hit my performance ... otherwise do I need to write the trigger for this...

What happens if the data increases in the table more than 10,00,000... because we are picking the last inserted data based on the device everytime?

Do i get any db schema for basic vehicle tracking?

Thanks in advance,
AJ


texts:
I'm trying to convert a CSV to PostgreSQL with ogr2ogr, but the CSV is coming from an API:

http://api.statbank.dk/v1/data/FT/CSV?delimiter=Semicolon&HOVEDDELE=*&Tid=*

I have tried following ogr2ogr:

ogr2ogr -f PostgreSQL PG:"dbname=postgres host=localhost port=5432 user=postgres password=xxx" CSV:"http://api.statbank.dk/v1/data/FT/CSV?delimiter=Semicolon&HOVEDDELE=*&Tid=*" -nln FT -lco SEPARATOR=SEMICOLON

But are getting

Unable to open datasource `CSV:http://api.statbank.dk/v1/data/FT/CSV?delimiter=Semicolon&HOVEDDELE=*&Tid=*' with the following drivers.....

I can get it working if I download the CSV and use ogr2ogr, but would prefer to fetch it straight from the API to reduce the steps in my script.


texts: Are there any ways to add a Shapefile layer using Geotools onto a map which created by Unfolding Maps? Or do I need to convert the shapefile(.shp) into other forms first?


texts: MySQL Polygon (Spatial) Insert Statement works sporadically

texts:
I know that MySQL doesn't have great spatial support at the moment, but is it that buggy that my insert statements would work sporadically? In general, longer coordinates seem to be failing on insert (there is no notice of them failing though, they just insert as NULL or show 0 rows affected on update). I believe though that I've seen large ones work and small ones fail. Here's my insert statement:

UPDATE places SET geom = PolygonFromText('POLYGON((-122.225255 37.890342,-122.22343 37.891445,-122.220923 37.89297...)),1') WHERE id=16

The only semi-logical things I can think of pertain to size/memory, or maybe invalid coordinates? Anyone have any ideas?


texts: After replacing table with updated data on Enterprise geodatabase, popups no longer work on Portal map

texts:
I have to do weekly updates on a table on SQL Server Enterprise geodatabase using Safe FME ETL Tool (drop/recreate). The table has the exact same name and exact same fields every time.
When I do replace the table, the features show up as needed, but the popup will not show on the portal map after the table update.
If I reshare the layers using ArcGIS Pro then the PopUps will show again, but this is adding another time-consuming step that I would like to eliminate.
I have tried multiple things, such as removing the layer and adding it again, also adding it "from web" using the rest url. Popups still do not show.


texts: Using AGSJSONRequestOperation in ArcGIS iOS Runtime SDK v100.2

texts:
In previous version of ArcGIS, the classAGSJSONRequestOperation had few methods to set target and actions. But in ArcGIS iOS runtime SDK v100.2 (latest version), there is no way to set the targets and actions.

Is it possible to use the AGSJSONRequestOperation class or any other replacements for this?


texts: Find a linestring closest to a given linestring

texts:
There must be already a solution to this problem - let's say we have five different linestrings. They might intersect or not. Given a sixth linestring, how do we find which one of the original five matches the sixth the closest, both in shape and space? A perfect match here would be if sixth linestring is exactly the same as one of the pre-existing ones. I the use case I'm looking at, one of the five would be almost exactly the same as the sixth one.

Anything in PostGIS, Java, Python would work for me, or even a general concept as well.


texts: Adding dynamic text to ArcGIS map documents

texts:
I've been making some map layouts and just noticed this command:

< dyn type="document" property="name"/>.mxd

It make possible to change automatically my current file name: http://prntscr.com/4sdzwl

I want to ask how does it works and how do I find some other commands that makes things like these quick and easy?


texts: ArcGIS 10.1 Archydrotools watersheed delineation AdjointCatchment error

texts:
The process of adjoint catchment processing considered during watershed delineation, it is complete but I only have a very small number of watersheds. I found the following error, in the report:

Generating adjoint catchments...
Error System.OutOfMemoryException: Exception of type 'System.OutOfMemoryException' was thrown.
   at System.Number.FormatInt32(Int32 value, String format, NumberFormatInfo info)
   at Microsoft.VisualBasic.CompilerServices.Conversions.ToString(Int32 Value)
   at ESRI.APWR.ApHydro.ApHydroConfig.MergeShape(ApProgworkMgr& pApProgMgr, IFeatureClass& pFCTracing, IFeatureClass& pFCSource, IFeatureClass& pFCTarget, ITable& pTableFS, String& NextDownFld, Int32& iMaxLoops, Boolean bMergeShape, Boolean bSaveRiverOrder, String FeatureIDFLD, ITrackCancel& trackCancel, IGPMessages& messages, Boolean bStopOnValue) in C:\Builds\HydroSolutions\10.3_ArcHydro\Sources\ArcHydroTools\src\ApHydro\ApHydroConfig.vb:line 4919
Pass 1: 416372 features processed, 387443 remaining. dt=394,3s dtTotal=423,0s
Error4 System.OutOfMemoryException: Exception of type 'System.OutOfMemoryException' was thrown.

I´m processing a large área (1.246.700 km²) and I have considered a stream threshold to a very small number 500. Is this a problem with processing capacity?


texts: Importing style information from .lyr file to Google Earth?

texts:
I have downloaded the geology data from the USGS website (http://pubs.usgs.gov/ds/781/OffshoreSantaBarbara/data_catalog_OffshoreSantaBarbara.html) and the .zip file comes with a .lyr file.

I have tried everything in my power to open this file on Google Earth - tried to convert it to .shp and .kml file using ArcToolbox - and nothing worked.

How can I use the style information contained in the .lyr file to visualize it on GE?

Edit: The problem I get when trying to convert to .kml is that I cannot find the saved file anywhere. It says that it was created, but I simply cannot find it.

When I first open the layer file it shows the red exclamation mark, but then I just add the .shp file and the image appears with the correct colors and symbols.

Then I try to convert using the ArcToolbox:

And I cannot find the file anywhere afterwards, although it says everything worked.