texts: Did the syntax for proj4 change with QGIS 3.10?

texts:
In QGIS 2.8 I was able to successfully do a transformation using the following proj4 syntax :

+proj=utm +zone=18 +ellps=clrk66 +datum=WGS84 +towgs84=-11.0584,163.8703,190.541,0,0,0,0 +units=m +no_defs +wktext

With QGIS 3.10 this no longer seem to work in the customised projection window.  Reading some of the documentation it seems QGIS 3.10 is now using proj6.  Is there a change that needs to be done in the syntax?


texts: I create a layout in QGIS, I'm going to put a grid (LAT and LONG); so the problem is my map scale is too small so the interval must be less than 1 degree.


texts: String quoting issues with sql using ST_PointFromText in psycopg2

texts:
I've been banging my head against this for a while, and am seeking some insight into what I'm doing wrong. I have the following snippet of code:

params = [PlaneRep.hex, PlaneRep.squawk, PlaneRep.flight, PlaneRep.isMetric,
    PlaneRep.mlat, PlaneRep.altitude, PlaneRep.speed, PlaneRep.vert_rate,
    dbdate, dbtime, PlaneRep.track, coordinates, PlaneRep.messages]
sql = '''
    INSERT into planereports (hex, squawk, flight, "isMetric", "isMLAT", altitude, speed, vert_rate, date_of_report, time_of_report, bearing, report_location, messages_sent)
    VALUES (%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, ST_PointFromText(%s, 4326), %s);'''
print(cur.mogrify(sql, params))
cur.execute(sql,params)

which yields the following error:

b'\n\t\tINSERT into planereports (hex, squawk, flight, "isMetric",     "isMLAT", altitude, speed, vert_rate, date_of_report, time_of_report,     bearing, report_location, messages_sent)\n\t\tVALUES (\'7c6d33\', \'3241\',     \'VOZ850  \', true, false, 9563, 822, 546.2016, \'2015-10-21\', \'04:32:41\', 235, ST_PointFromText(\'POINT(150.464266 -34.754929)\', 4326), 133);'
Traceback (most recent call last)
File "./jsonfiddling.py", line 92, in <module>
logToDB(dbconn, plane)
File "./jsonfiddling.py", line 34, in logToDB
 cur.execute(sql,params)
 psycopg2.ProgrammingError: function st_pointfromtext(unknown, integer)    does not exist
LINE 3: ...63, 822, 546.2016, '2015-10-21', '04:32:41', 235, ST_PointFr...
                                                         ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

Any ideas of what I'm doing wrong?


texts: Raster slightly displaced after exporting from PostGIS database

texts:
With the following queries I've exported a raster (GTiff) from my postgis database:

1.Creating a large object:

WITH buffer AS (
    SELECT
        id,
        ST_Buffer(geom,50) AS geom
    FROM polygon_table),

clip AS (
    SELECT
        partnumber,
        ST_Union(ST_Clip(a.rast, 1,b.geom, TRUE)) AS rast
FROM aspect a
        JOIN buffer b
        ON ST_Intersects(a.rast, b.geom)
            GROUP BY partnumber)

SELECT
    oid,
    lowrite(lo_open(oid, 131072), png) As num_bytes
    FROM (
        VALUES (lo_create(0),
               ST_AsGDALRaster(((SELECT ST_Union(rast) FROM clip)),'GTiff')
           )) As v(oid,png)

2.Exporting it to a GTiff file (e.g. with oid=176148)

SELECT lo_export(176148, '/tmp/demo_rast.tif')

(3. Deleting the large object: SELECT lo_unlink(176148))

The aspect raster data in the postgis database has the following information:

 rid | scalex | scaley |  srid
-----+--------+--------+-------
  1  |   25   |  -25   | 31468

The aspect raster resolution is 25m. Tiling is 100x100.

When load the raster into QGIS it is slightly displaced. The colored aspect raster is the original. The grayish one is the exported one.

The amount of the displacement is about 2.65m to the west in x direction and about 1.35m to the north in y direction.

Where does this displacement come from?

UPDATE

The problem appears exclusively in QGIS.
When I load both rasters into GRASS GIS there is no displacement. I didn't used the option Override projection (use location's projection) when I import the rasters into GRASS GIS.
As you can see here:

So is there a general problem (bug) which appears here in QGIS???


texts: Drawing the most efficient path between multiple points

texts:
I am trying to find the most efficient path between multiple points on a map.
A bit of background to help understand what I need. I am making a printed map that has LEDs poked through from behind and each LED has a wire 10cm between the next. I am trying to use as many LEDs as possible and minimise wasted LEDs left behind the map.
I have all the points as coordinates in QGIS.
Is there a way to have QGIS draw a line of maximum "10cm" or just a most efficient path between all the points?


texts: Estimate missing coordinates in a trajectory of points in Python

texts:
I have points with longitude, latitude and date, being the trajectory of a boat taken every day.
One point is missing and I have to estimate it with Python.

Is there built-in libraries doing that or should I only consider that:

boat is going at constant speed that I compute from other points
boat is keeping the same direction

and compute the point from that ?

--edit--

To have a good estimate I would like to be sure that the estimated point is not on the ground.


texts: Landsat images with more than 3 colored bands

texts:
I’m looking for a website where I can find Landsat images (or similar to it) with more than just 3 color bands. http://glovis.usgs.gov/ only gives me 3 bands. Is there a way to obtain more band?
Thanks for your answers!


texts: Visualization of satellite image - directly in Python

texts:
Is there any way (a package or what not), to visualize a GeoTIFF image directly in Python?

I don't want to use ENVI or ArcGIS or anything else.


texts: Using both ArcGIS API for JavaScript 3.16 and 4.0

texts:
I use ArcGIS API for JS 4.0 but I didn't find a widget to add a scalebar to my web application.
So, I tried to use ArcGIS API for JS 3.16 which has a scalebar widget, but I don't know if it's a good idea to mix these different versions.


texts: Unexplained coordinates mutation with readOGR in R

texts:
I am trying to intersect two spatial objects in R (point and polygon) while preserving the contained information in the dataframes of both.

However, the intersection produces an empty polygon, while both objects do perfectly overlay in QGIS 2.16.2.

After I checked the coordinates of the spatial polygons (test_polys), I find them strangely altered, such as coordinates like  -44524.15 , 167543.0.

My efforts to change the coordinate system (Step 3) did not succeed.
Some ideas how to fix it and explanations, where the coordinate modifications could originate from?

I use R.Studio Version 1.0.143 and rgdal 1.2-8 as well as raster 2.5-8.

Please find the input files in this dropbox folder

And the code I use:

# 1) Initialization -----------------------------------------------------

library(rgdal)
library(sp)
library(raster)

setwd("SET YOUR DIRECTORY")
list.files()

# 2) Read data ----------------------------------------------------------

test_polys  <- readOGR(dsn = "C:/YOUR PATH/Desktop/Stackexchange", layer = "Tests_10subsections")
test_points <- readOGR(dsn = "C:/YOUR PATH/Desktop/Stackexchange", layer = "calc_centroids")

test_polys@polygons # Strange coordinates already at this step, such as "-44524.15 167543.0"

# 3) CRS Homogenization -------------------------------------------------

proj4string(test_polys)       <- CRS("+proj=longlat +datum=WGS84 +no_defs  +ellps=WGS84 +towgs84=0,0,0")
# ERROR "Error in `proj4string<-`(`*tmp*`, value = <S4 object of class "CRS">) : Geographical CRS given to non-conformant data: -45033.0103 167554.5511"

test_polys@proj4string        <- CRS("+proj=longlat +datum=WGS84  +no_defs +ellps=WGS84 +towgs84=0,0,0")
test_points@proj4string       <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

    test_points  <- spTransform(test_points, CRS(test_polys))
# Error in CRS(test_polys) :
# no method for coercing this S4 class to a vector
# In addition: Warning message:
# In is.na(projargs) : is.na() applied to non-(list or vector) of type   #'S4'

# 4) Intersection while preserving point data adding polygon data -------

intersected  <- raster::intersect(test_points, test_polys)
# EMPTY "data" with: [1] id         Roof_area  DTMN11     FR11       SEC11      SS11       BGRI11     LUG11      LUG11DESIG
#<0 rows> (or 0-length row.names)


texts: How to add Tiles made from an OSM PostGIS Database to QGIS

texts:
Following this tutorial I was able to configure a tile server correctly using a PostGreSQL/PostGIS database, osm_tiles, renderd/mapnik and osm2pgsql in CentOS 6.

The goal of this is to be able to open tiles that have been uploaded to my DB in QGIS and some interfaces like Leaflet/OpenLayers.

Link between tile server and QGIS  enables me to visualize the "World tile" that shows the server is working, using the URL x.x.x.x/{z}/{x}/{y}.png in QGIS.

Now I would like to be able to "read" the tiles that are in my PostGIS DB (I uploaded the Benin) but have no idea how to do it...


texts:
I have downloaded VIIRS day/night Band data from here but when I try to open tiff files in QGIS software, it doesn't show like usual night light files and it is similar to cloudy sky image (see the screen shot below).

What should I do to open and see it truly ?


texts: QGIS--computing local standard deviation (from many bands) to each cell

texts:
I have HDF images with 365 bands each. Is there a way to compute local standard deviation across all of these bands with results appearing in each cell?

I have tried Zonal Statistics (with polygons converted from each raster cell as the input vector layer) but it does not appear to work for multiple bands. I also tried using SAGA via QGIS but it gave an error stating multi-band functionality was not supported.

Any one know of a way I can compute this using Raster Calculator or GRASS's r.series?


texts:
I am getting warm with post PostGIS.
I am working in the Command Line. What I am wondering is how do I "get" the results of PostGIS.

So imagine I have this

select name, neighname from neighborhoods where ST_Intersects(geom, ST_GeomFromText('LINESTRING(586781.701577724 4504202.15314339,586863.51964484 4504215.9881701)', 26918));

That will prompt me something like this in the command line

    name    | neighname
------------+----------
   Mainst   | West Vill
         (1 row)

So far, so good. But how would I proceed working with this selection?

Can I somehow save that as a new shape?

Can I get the results into QGIS without the step of saving it as a shape (if that is possible)?

How can I proceed working with that selection in PostGIS?


texts: Importing/copying geodatabase schema from multiple tables to single table?

texts:
I have an existing ArcGIS 10 geodatabase which has a complicated schema containing multiple tables and joins, and I wish to simplify this for faster web delivery.

For example, Featureclass1 is related to Table1, Table2, TableN, and I wish to copy or import the attributes from each Table onto Featureclass1. This will result in a single Featureclass containing all attributes.

What is the fastest way to do this in a repeatable manner (ie, without doing it manually in ArcCatalog)?

This is what I'm considering:

When creating a new featureclass or table in ArcCatalog, there is the option to Import the schema from an existing table/FC. However, this only allows you to import from Table1 (if you then choose Table2, it drops the fields from Table1).
The Geodatabase Diagrammer is a great tool for reporting on an existing schema, and designing a new one - but I can't see any options to import multiple fields simply (it seems to require manual entry of each new field).
writing a Python script to do it iteratively. Is this the best option?

Note that I'm only concerned with the schema at this stage, not the actual data import.


texts: Combining corner-touching polygons into a single feature

texts:
I have created a polygon layer from rasters, and only have a single classification ID (and calculated acres) for each polygon.  I would like to "dissolve" all touching features into larger, individual polygons.  The 3 areas circled in red (below) contain green polygons/pixels that only touch on the points of the squares.  Intuitively (from an end user perspective) all of the green polygons within a circled area should be connected and treated as a single polygon because there are no "breaks" where the features/cells do not touch.

What would be the best method to accomplish "joining" or "dissolving" all of these features using ArcGIS.  I have ArcInfo and Spatial Analyst extension.  If there is a better methodology for doing this in Spatial Analyst, I can pre-process the polygons that way, too.


texts: Getting inconsistent raster output from gdal_calc.py raster calculation

texts:
I'm attempting to do what is generally a straightforward raster calculation when I do this in QGIS, but am trying to do this with gdal_calc.py in terminal. The output is close, but produces banding and inconsistencies in parts of the output raster. I'm wondering if there is anything I'm overlooking. Code to calculate the image is below:

gdal_calc.py -A <input>  --outfile = <output> --calc="(A/nanmax(A)) * 65535" --NoDataValue=0


texts: Saving edited vectors to new file does not save changes?

texts:
Perhaps I do not understand how QGIS saves edits but here is what I am experiencing.

Load up a vector .shp file
Immediately "SAVE AS" and give it a different name.
Delete some vector lines on the new file.
Save the edits.
Drop the edited .shp file into Global Mapper.  All the deleted vector lines are still there!

What am I doing wrong or misunderstanding on how QGIS saves edits?


texts:
I am keen to make a career in GIS.

What should I work upon to possess the adequate skills and can you recommend some good colleges for a GIS-oriented Postgraduate course?


texts: Scaling factor for MODIS (MY09GA 006 or MYDTBGA 006) product for timeseries data

texts:
I was using Modis data (MYDO9GA) product for my research and I was following this document as a user guide.

https://landweb.modaps.eosdis.nasa.gov/QA_WWW/forPage/user_guide/MOD09_UserGuide_v1.4.pdf

I applied the scaling factor (0.0001) and valid range filter (-100 to 16000  from document) on my albedo data (from MYDTGA006 product) but  values do not seem normal. For example  I am getting values between -0.12 to 0). However albedo range should be between 0-1. I was not using raster data. I extracted time series data based on lat long using following link.   https://lpdaacsvc.cr.usgs.gov/appeears/

I think scaling is not necessary for timeseries point data.For example solar zenith data looks fine without scaling factor. I was wondering if anyone has experience with his data product timeseries.

Also link to the question asked previously to provide context. Again I am not using raster image (like used in the question below). I am using timeseries data.
Reprojected MODIS NDVI has range from -32768 to 32767, expected -1 to 1


texts: get "Cannot find coordinate operations" in GDAL 3; works in GDAL 2

texts:
I get the following error when trying to use GDAL 3.1.2 for computing the LLA from my custom gnomonic projections. Am I setting it up wrong? Missing some information/files? Would it be a misconfigured build? It worked fine for me in GDAL2.
Cannot find coordinate operations
from `PROJCRS["Customized Gnomonic",BASEGEOGCRS["Customized WGS84",DATUM["Custom Datum",ELLIPSOID["WGS84 Radius + Elevation",6379509.75,inf,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433,ID["EPSG",9122]]]],CONVERSION["Gnomonic",METHOD["Gnomonic"],PARAMETER["Latitude of natural origin",41.7480066139774,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",-112.023800209591,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
to `EPSG:4326'

Basic code to set it up:
var sr = new SpatialReference(null);
sr.SetProjCS("Customized Gnomonic");
sr.SetGeogCS("Customized WGS84", "Custom Datum", "WGS84 Radius + Elevation",
            _earthRadius + origin.Altitude, double.PositiveInfinity, null, 0.0, null, 1.0);
sr.SetGnomonic(origin.Latitude, origin.Longitude, 0.0, 0.0);


texts: ArcPy inconsistently accepts SQL queries for MakeFeatureLayer or SelectLayerByAttribute

texts:
I need to buffer a subset of features in a particular feature class. To do this, I am first creating a feature layer (using a SQL query to select only what I want), and then running Buffer. So my ArcPy script includes this little function to buffer my handful of feature classes:

def funcBuffer(inFeature, outFeature, bufDist, querySQL):
    arcpy.MakeFeatureLayer_management(inFeature, "tempLayer", querySQL)
    arcpy.Buffer_analysis("tempLayer", outFeature, bufDist, "FULL", "ROUND", "ALL", "")

This works perfectly on most of my feature classes. However, there is one that fails. The Python error is:
ERROR 000230: Failed selecting with "OH_UG" = 'UG'

The help page for that error suggests that the SQL query is invalid. However,

the field OH_UG exists, and there are some features with the value "UG"
it successfully validates in Query Builder
I can use it to Select by Attribute with no problem (python snippet from within ArcMap: arcpy.SelectLayerByAttribute_management("TransmissionLines","NEW_SELECTION","OH_UG = 'UG'"))
similarly structured queries work on other feature classes

I've tried different ways of passing the SQL query and all fail: '"OH_UH" = \'UG\'', 'OH_UH = \'UG\'', """ "OH_UG" = 'UG' """

I revised the script to do the selection after making the feature layer:

def funcBuffer(inFeature, outFeature, bufDist, querySQL):
    arcpy.MakeFeatureLayer_management(inFeature, "tempLayer")
    arcpy.SelectLayerByAttribute_management("tempLayer", "NEW_SELECTION", querySQL)
    arcpy.Buffer_analysis("tempLayer", outFeature, bufDist, "FULL", "ROUND", "ALL", "")

and now MakeFeatureLayer works, but SelectLayerByAttribute throws ERROR 000358: Invalid expression -- still a complaint regarding invalid SQL.

I attempted to select using a different (numeric) field in the table, and am still getting the same 000358 error.

So I am not sure what's going wrong -- is it the SQL, is it the particular feature class, is it gremlins? How can I work around the problem in ArcPy?


texts: Band Names in Google Earth Engine - Modis LST data

texts:
I am working with some modis LST data, and am running into the following problem:
I am running two functions on the data, which I will need to combine later. For that reason, the band names need to be the same. When I run the code below however, one output has band [LST_DAY_1km], and one has band [constant]. This obviously doesnt add up, and I am confused where the problem is coming from, given I'm using the same input data.
See code below:
    // Use the expression function to generate scaled LST.

var scaledlst1 = lst.map(function (lst) {
  return lst.expression(
    '(lst/30)',{
      'lst': lst.select('LST_Day_1km')
    });

});

Map.addLayer(scaledlst1, {band:'LST_Day_1km' ,max: 1 , min: 0}, 'scaledlst1');

//

var scaledlst2 = lst.map(function (lst) {
  return lst.expression(
    '(2.5-(0.05*lst))',{
      'lst': lst.select('LST_Day_1km')
    });

});

Map.addLayer(scaledlst2, {band:'LST_Day_1km', max: 1 , min: 0}, 'scaledlst2');

Does anyone have any suggestions as to why this might be?
My next step is to take the minimum value per pixel between the two products:
//select images

var ListOfImages1 = scaledlst1.toList(scaledlst1.size())
var ListOfImages2 = scaledlst2.toList(scaledlst2.size())

var LST = ee.Image(ListOfImages1.get(0))
var LST2 = ee.Image(ListOfImages2.get(0))

//Take min value
var scaledLST = ee.ImageCollection.fromImages([LST,LST2]).min()

Map.addLayer(scaledLST, {max: 1 , min: 0}, 'Scaled LST');

Which yields the following error:
Scaled LST: Tile error: Expected a homogeneous image collection, but an image with incompatible bands was encountered:
First image type: 1 bands ([LST_Day_1km]).
Current image type: 1 bands ([constant]).
Image ID: 1
Some bands might require explicit casts.


texts: Set a WMS URL in QGIS using project variables

texts:
I would to know if it is possible to set a WMS url in which I would use project variables as parameters?

Let's say I have:
QGIS project variable : my_variable = some_query
WMS url : https://wxs.myserver.com/wms?LAYERS=OI.OrthoimageCoverage& EXCEPTIONS=text/xml&FORMAT=image/png&SERVICE=WMS&VERSION=1.3.0& REQUEST=GetMap&STYLES=&CRS=EPSG:4326&WIDTH=256&HEIGHT=256

I want to replace the LAYERS=OI.OrthoimageCoverage with my predefine project variable as follow:

https://wxs.myserver.com/wms?LAYERS=@my_variable& EXCEPTIONS=text/xml&FORMAT=image/png&SERVICE=WMS&VERSION=1.3.0& REQUEST=GetMap&STYLES=&CRS=EPSG:4326&WIDTH=256&HEIGHT=256

Is there a way to do this?


texts:
I need to access to the function ee.classifier.cart to understand all the passages of the codes and possibly to modify it for my purposes.

Is there any way or I need to be an Earth Engine Developer?


texts: How can I strip out Z values in PostGIS?

texts:
How do I get rid of the Z dimension in PostGIS?  Specifically, I have some MultiPolygonZs that I want to copy into another table as MultiPolygons.  I don't care about the Z values - it's okay to just throw them away.

There must be a simple way, but I've been unable to find any functions that do this...  Any tips are greatly appreciated!


texts: QGIS to Unreal Engine (Possibly via Terraform Pro [Free]) Totally Stuck!

texts:
Totally in a tiz & getting really confused (as not a GIS sort), but I am trying to ideally get Full Topography data into Unreal Engine using QGIS, Ideally per whats in the OpenTopoMap & Ocean, depth/detail say like in Esri Ocean/sailing/marine chart & inland water depth/streams/rivers. Or better as did try downloading DEM Height data to then get even more confused.

This is likely what I will try to use once I can figure how to process all I need in QGIS, into Unreal Engine, I also have L3TD if any use, the tutorial, which will likely make lots of sense to you clever people....
 https://terraformpro.com/tutorials/tutorial-1/introduction/

As its an open source, free educational project/interactive story idea, even better would be someone volunteering to help with the know-how and who can output all the layers I´d need for 1M scale (if possible for the area Around & including St Martin in the Carribean. N018W064, N018W063, N017W064, N016W063.)

Obviously, if anyone could help directly I will credit with heaps of praise and link to our sites, social media as I will in regards to the sources of the Data used in the files created.

I myself am Autistic and have ADHD and the project idea is for kids & adults like me. The idea is to have a physical book, story that links into the interactive story/game. The interactive story in book and UE4 is linked to the maps in the covers of the book, front Ecology, back Weather more storms etc as set a little bit in the future where global warming has really hit. But dealing with the Geography/Ecology bit first, second weather).


texts: How can I manipulate the timestamps in a GPX file?

texts:
I have a GPS tracked stored in a GPX file. Now I'm looking for a tool that allows me to set the start time of the track and the duration, and the tool should automatically evenly distribute the appropriate timestamps for all my recorded points in the track.

The tool should be free and run on OSX.


texts: I've downloaded the bonn version of QGIS. Using the 'regular points' vector analysis tool, I've set some 5m regular points over a field so that I can use the x,y coordinates to plot them and measure their z coordinate in the field with a gnss. I am gathering height data for a DTM. When I export the layer in QGIS as a csv file it doesn't give me the coordinate x,y values, only the id.


texts: Properly deal with broken project data source Analyze Errors when publishing Geoprocessing Service (Python / ArcGIS Server 10.2)

texts:
When publishing a geoprocessing service, I noticed that arcgis does not like to see:

os.path.join(os.path.dirname(__file__), 'projections')

This line seems to throw an broken project data source while analyzing the service for publishing to arcgis server.

I'm guessing its because the only the main .py file will be copied into the services directory and Esri is giving me a heads up that the data source is going to break because 'projections' folder will NOT be copied over.

In my case, I also have some unit test which reference a test_data folder and shouldn't affect the tool when running off the server. Unfortunately I still can't publish.

To publish, I:

Moved my relative path references to a separate module named service_config.py
imported the path constants from service_config.py
Added a .pth file to Lib/site-packages to allow the services directory version of the script to still have references to the projection folder without needing it to be copied to the arcgisinput directory.

I feel like my solution is a hack. Any suggests on how to improve this tool configuration?

My code can be found here:

https://github.com/brendancol/esri-extract-data-where


texts: SQL date Query of personal geodatabase in ModelBuilder with inline variable

texts:
In ModelBuilder I'm trying to create a table from a subset of dates in a Personal geodatabase  table. I'm using "Make Table view" tool.

The date field in the table is [DateCreated]. I have two model parameter variables %Start Date% and %End Date%, both with data type set as date. I would like to return values from the table between these two dates.

I'm struggling to find resources to describe the required syntax. I'm also not very familiar with SQL.

So far I've tried the below with no success

[DateCreated] >  #%Start Date%# AND [DateCreated] < #%End Date%#

[DateCreated] >  %Start Date% AND [DateCreated] < %End Date%

[DateCreated] >  '%Start Date%' AND [DateCreated] < '%End Date%'

In all of these example the tool returns all rows


texts: Adobe Flex layer is not displayed, but can zoom to it

texts:
I am using Adobe Flex to develop a widget to allow the user to turn on/off 2 layers separately. The strange thing is only 1 layer can be displayed on the map when the layer is turned on, the other layer does not appear when turned on (but can still zoom to the layer). I cannot figure out what happens. The difference of the 2 layers is that the layer not appearing is joined to a database table. Will that affect it?

I am using Adobe Flash Builder 4.6 and ArcGIS API for Flex 3.6.

Any ideas on it?


texts: Polygonize RGB mage considering all 3 bands simultaneously

texts:
I'm using the rasterio.features.shapes method to generate GeoJSON information from a PNG file. This answer shows what I'm trying to do, except I need to analyze the RGB bands simultaneously otherwise the polygons won't be correct. I'm circumventing this by converting first the triplet of 3 int8 to a single 32bit integer using the bitshift operator, but I feel I'm missing the obvious solution to combine the 3 bands together in some way. I'm also using PIL and numpy to process the image.
import numpy as np
from PIL import Image
from rasterio.features import shapes

def find_polygons():
    im = Image.open(f"myimage.png")
    arr = np.array(im)
    out = np.vectorize(rgb_to_int32, otypes=[np.int32])(*np.rollaxis(arr, 2, 0))
    for i, (s, v) in enumerate(shapes(out)):
        r, g, b = int32_to_rgb(v)
    print(i)

def rgb_to_int32(r, g, b):
    return (r << 16) + (g << 8) + b

def int32_to_rgb(n):
    n = int(n)
    return n >> 16 & 0xFF, n >> 8 & 0xFF, n & 0xFF

This produces the desired result, but I feel like I'm missing the obvious solution. Any tips in the right direction?
For reference, this is the image I'm trying to extract data from:


texts: iterating and exporting all combinations of features as shapefiles

texts:
I am trying to export individual shapefiles of all 2 pair combinations from a multi-feature polygon shapefile.

To explain this better as set of FID's 0, 1, 2, 3, 4 would be exported as individual shapefiles of FID pairs:

('0', '1'), ('0', '2'), ('0', '3'), ('0', '4'), ('1', '2'), ('1', '3'), ('1', '4'), ('2', '3'), ('2', '4'), ('3', '4')

I am not so good with python so my attempts with itertools.combination etc have not worked out.


texts: QGIS How to find QVariant of existing attribute

texts:
I want to add attributes, without features, that exists in one layer but not in the other using the code below. I see that I am not obtaining the "QVariant" from anywhere and just assumes it will be Double. Is there a way to obtain the correct information from the source layer? Or is there a better way of doing this?
diff = list(set(ly1.fields().names()).difference(ly2.fields().names()))
for att in diff:
    res = ly2.dataProvider().addAttributes([QgsField(att, QVariant.Double)])
    ly2.updateFields()

I have looked into using:
type = ly1.fields().field(item).typeName()

This however returns, e.g., int8. But unsure how to use it.


texts: ERROR 001269 : Compressing the service definition failed.ArcGIS 10.1,Python 2.7

texts:
I am trying to publish an mxd using python script. I invoke this pyhton script from c# application, after creating the sddraft file, the system throws up this error "001269 : Compressing the service definition failed."
I looked up online and it says that it gives this kind of error with 64 bit python, I was using python for visual studio(not Iron Python). And it also says online that "ArcGIS for Desktop was still using the 32 bits"
Is there a way to publish the mxd or anything else which I should consider? Any help on this is most welcome.


texts:
I have kind of the same question as "greyline" : How to implement 3NF (third normal form) in ArcGIS attribute tables?

I need, among other things, to place foreign key in my tables to combine them and I don't know how to do it with ArcGIS.
I tried at first with "New -> Relationship Class" but I didn't have control on my attributes to select properly and avoid redundant data (I don't know the equivalent of "distinct" in ArcGIS commands).
I am now trying to create a model that could combine the tables, but I don't know how to specify that each row for my "Date-t" tab, has to be linked to a single proper table in my database?

My model's structure:

"Make Table View" details:

"Add join" details:

I am ready to change all the method if mine is wrong. I just cannot use the PostGIS since I don't have it.

The aim of this operation is to be able to fetch data from my tables using Python, but I think it's an other question.


texts: Post Orthorectification error reduction?

texts:
I am trying to orthorectify an old aerial photo that has control points. When I run an output and then validate it I find some offset or distortion.

What should I do to fix this?

Does anyone have an example of pre and post orthorectification to help me?

I also seek any additional sources of information on this process that anyone can share.


texts: How to make an instance of a GUI every time when click button

texts:
I create a custom plugin in QGIS 2 with Qt Designer. When I press button of the plugin I want to open a dialog with a table that has loaded data from the db. The problem is that I am confused on how to call GUI class to make an instance of the GUI every time that I click the plugin.
I have this class that creates the GUI, but I am confused on how to use it.
class myClassGui(self):
    def setupUi(self, Dialog):
        Dialog.setObjectName(_fromUtf8(myGui))
        Dialog.resize(911, 540)
        self.btn= QtGui.QDialogButtonBox(Dialog)
        self.btn.setGeometry(QtCore.QRect(510, 470, 341, 32))
        self.btn.setOrientation(QtCore.Qt.Horizontal)
        self.btn.setStandardButtons(QtGui.QDialogButtonBox.Cancel|QtGui.QDialogButtonBox.Ok)
        self.btn.setObjectName(_fromUtf8("btn"))
        self.myTable= QtGui.QTableWidget(Dialog)
        self.myTable.setGeometry(QtCore.QRect(30, 20, 841, 411))
        self.myTable.setRowCount(1)
        self.myTable.setColumnCount(18)
        self.myTable.setObjectName(_fromUtf8("myTable"))
        item = QtGui.QTableWidgetItem()
        self.myTable.setHorizontalHeaderItem(0, item)
        item = QtGui.QTableWidgetItem()
        self.myTable.setHorizontalHeaderItem(1, item)
        item = QtGui.QTableWidgetItem()
        self.myTable.setHorizontalHeaderItem(2, item)
        item = QtGui.QTableWidgetItem()
        self.myTable.setHorizontalHeaderItem(3, item)
        item = QtGui.QTableWidgetItem()
        self.myTable.setHorizontalHeaderItem(4, item)
        self.retranslateUi(Dialog)
        QtCore.QObject.connect(self.buttonBox, QtCore.SIGNAL(_fromUtf8("accepted()")), Dialog.accept)
        QtCore.QObject.connect(self.buttonBox, QtCore.SIGNAL(_fromUtf8("rejected()")), Dialog.reject)
        QtCore.QMetaObject.connectSlotsByName(Dialog)


texts:
I am getting some weird errors while trying to implement multithreading with ArcObjects.
I am trying following things on per thread basis:

1) Open SDE Workspace
2) Open Feature class
3) Perform Icursor

Eventually, It fails with some errors like:
"Memory could not be read or write" or "AccessViolationException".
Looks like there are some memory leaks (though I don't know exactly where).

For more details please refer to the code below:

foreach (string featureClassName in settings.FeatureClassNames)
{
            Thread thread = new Thread(new ParameterizedThreadStart(ReadFeatureClass));
            thread.ApartmentState = ApartmentState.STA;
            object[] parameters = CreateParameters(featureClassName);
            thread.Start(parameters);
}
    private void ReadFeatureClass(object parameters)
    {
        Settings settings;
        string featureClassName;
        ExtractParameters(parameters as object[], out settings, out featureClassName);
        IFeatureWorkspace featureWorkspace = OpenFeatureWorkspace(settings);
        ITable featureClass = featureWorkspace.OpenTable(featureClassName);
        Marshal.ReleaseComObject(featureClass);
        Marshal.ReleaseComObject(featureWorkspace);
      }

    private IFeatureWorkspace OpenFeatureWorkspace(Settings settings)
    {
        IWorkspaceFactory workspaceFactory = new SdeWorkspaceFactory();
        IPropertySet connectionProperties = new PropertySet();
        connectionProperties.SetProperty("SERVER", settings.Server);
        connectionProperties.SetProperty("INSTANCE", settings.Instance);
        connectionProperties.SetProperty("USER", settings.User);
        connectionProperties.SetProperty("PASSWORD", settings.Password);
        connectionProperties.SetProperty("VERSION", settings.Version);

        IFeatureWorkspace featureWorkspace = workspaceFactory.Open(connectionProperties, 0) as IFeatureWorkspace;
        return featureWorkspace;
    }


texts:
I am interested in creating a list of all the IDs from Table_A and whether or not they intersect with any of the geometries from Table_B.

The following would work:

SELECT Table_A.id, bool_and(st_intersects(Table_A.geom, Table_B.geom)
FROM Table_A, Table_B
GROUP BY Table_A.id

The problem is that this means checking every geometry in Table_A against every geometry in Table_B which is very intensive. How do you query until you find the first intersect, return TRUE and then move onto the next geometry in Table_A?

The extension to this, is to then check each row in Table_A for an intersect in several other tables, with the result as follows:

Table_A_ID | Table_B_Intersects | Table_C_Intersects | ...
0001       | True               | False              | ...
...        |  ...               | ...                | ...

This answer How to find all features in one table that intersect features in another table in PostGIS? from comments helps with the single case, but I cannot see how to apply it in this instance.

Edit: moved from comment, tried the following, but crashes pgAdmin:

 SELECT Table_A.id,
 (CASE WHEN EXISTS
    (SELECT st_intersects(Table_A.geom, Table_B.geom))
 THEN TRUE ELSE FALSE END) AS Table_B_Intersects,
 (CASE WHEN EXISTS
    (SELECT st_intersects(Table_A.geom, Table_C.geom))
 THEN TRUE ELSE FALSE END) AS Table_C_Intersects,
 ...
 FROM Table_A, Table_B, Table_C, ...


texts: I have a 1-m continuous DEM and a vector grid with ~160,000 250m cells. For each cell in the grid overlaying my DEM, I want to retrieve the 50th percentile pixel value. I'm curious if there's a way to retrieve this percentile value using zonal stats in R?


texts: How to create polygons with straight lines and right angles in QGIS?

texts:
When creating a polygon, how do I do it so that the lines are straight (90, 0 degrees, etc.) and when I create corners they can be right angles?
For instance, how do I create a perfect square or rectangle?


texts:
I am trying to use right to left language in CesiumJS label, but the result is reverse ordered. What can I do?

 viewer.entities.add({
        position : Cesium.Cartesian3.fromDegrees(-75.59777, 40.03883),
        point : {
            pixelSize : 10,
            color : Cesium.Color.YELLOW
        },
        label : {text : "מה המצב" }
    });

Yields image:


texts: How to get the most current data from a Temporal Feature Service?

texts:
I have created a Temporal Feature Layer on our ArcGIS Server, which records all the locations, along with timestamp & other attributes of moving vehicles.

I used the TimeSlider Sample, to create a custom App which shows this data along with a TimeSlider.

When user presses play, it shows all the locations for a Given Vehicle (in that time extent) on the Map.

Is there a way to configure the timeSlider, or even make a custom query to the REST Endpoint of the FeatureService, to get the most current data/location for each vehicle in that particular time extent, and not all?

(I'm basically trying to get this kind of behavior in a web app)


texts: Layer not loading into ArcGIS Pro Catalog window

texts:
Using ArcGIS Pro 2.7, Learn ArcGIS trial license
I imported a feature class vector layer created in ArcMap into ArcGIS Pro. Under layer properties, I was able to copy the data sources metadata to the layer and access the metadata description in the Catalog window. I did this by navigating to the project map in the Catalog window, and clicking on the layer name listed on the left side of the window.
As part of an existing project in ArcGIS Pro, I re-projected the vector layer as a new feature class in the map project, however the new layer is not appearing in the layer list in the Catalog view. At the bottom of the Catalog view the status is Loading..., apparently unable to upload the new layer into the Catalog view even though the new layer exists in the contents pane of the map project.


texts: How to truncate extension off in-line variable %name% in ModelBuilder?

texts:
I have seen similar questions but no clear answer and I am afraid that this can't be done in ModelBuilder alone.

Is there a way to truncate the extension off the %name% variable? I am working with .tif and .asc files that I use to resample and do zonal statistics but since they have file extensions the tools throw an error when using %name% and I am forced to use %n% and then rename the files using R. I have tried "%name%"[x:y] as suggested in another question but doesn't work.

Any suggestions? This is ArcMap 10.2 but had the same problem in 10.0


texts:
How can I move a polygon to specific place. Vertex should be join with one another.
here is a picture


texts: Air corridor in spherical surface: calculating great circle tangent points with small circles

texts:
I have the necessity to draw on a sphere (earth) an air corridor.

In the cartesian plane, this corridor is made by some circles and from outer tangents that connect these circles as you can see in this image (central part of the image):

In general circles can have different radii. In flat Cartesian plane I can calculate the two outer tangents (when possible) of two circles, using for example formulas that I can found in this page.

My problem is a bit more complicated. I must to perform same calculus in spherical surface. You can see more clearly in this image:

In this case I've two small circles and I want to calculate tangents between these two circles in spherical surface.

Onto a sphere I don't have straight lines, but geodesics, and in this case they are represented by great circles. I am able to calculate a great circle trajectory between two points using the Aviation formulary.

What I miss is a way to find the small circles tangent points, in order to have the initial and starting points for my geodesic tangents.

How can I calculate them?


texts: Why does adding polygon in QGIS sometimes create a feature with edge too far across the screen? What am I doing wrong?

texts: While digiting/creating new polygons in QGIS some polygons have elongated edges. Some features have no geometry (area=zero) too. What could be the reason?  I use 3 Monitors, one for QGIS, one for GoogleEarth and one of Laptop. The version of QGIS is 2.18.11


texts: Select points located within Multi Ring Buffer

texts:
I have two point datasets, 1 showing the location of road accidents, the other showing the location of traffic intersections. I have created a Multi Ring Buffer at 25m intervals around the traffic intersections buffer (ie: a buffer at 25m, 50, 75m, etc.)

I want to select points from the road accidents buffer that fall within the 25m buffer, excluding those outside.

I then want to select points within the 50m buffer, excluding those which fall inside the 25m buffer and outside the 50m buffer. I'll then repeat this process for the remaining buffer.

How would I go about selecting points within the 50m buffer while excluding those that fall outside it and within the 25m buffer?


texts: "Transferring" ArcGIS Online content from one organization to another

texts:
My organization is currently brainstorming with another organization on how continue providing access to about a half of dozen web applications and data through ArcGIS Online.

Currently, the organization we are in talks with wants to be able to hand off the content/web applications on ArcGIS Online to us to host on our ArcGIS online organizational account.

I haven't found anything that steers us in the right direction on GIS SE and on ArcGIS Online FAQ here appears that this feature isn't entirely possible.

How would we go about doing this? If at all possible?


texts: Warp MODIS Blue Marble Projection to EPSG:102009 QGIS

texts:
I have been trying to reproject the land_shallow_topo_west.tif file for MODIS from the following site: https://visibleearth.nasa.gov/images/57752/blue-marble-land-surface-shallow-water-and-shaded-topography
But I can't seem to get it to work in QGIS, I keep getting the error:
ERROR 1: The transformation is already "north up" or a transformation between pixel/line and georeferenced coordinates cannot be computed for land_shallow_topo_west.tif. There is no affine transformation and no GCPs.
When I try to warp the TIF to EPSG:102009. I am not sure if I have to assign some other projection first to the image, I am assuming the default is WGS84? However none of the reproductions I've tried have been successful. Does anyone know the steps in QGIS to get this to work?


texts: FME - Possible to return a value from workspace?

texts: I'm calling my workspace from the command line (either via Python or .NET to a system/os call). I'm parsing DWGs and pushing attributes to SQL Server in the workspace. Either the data validates in the workspace or it doesn't. I need to return a value of 'pass/fail' to my calling application from my workspace based on whether or not the data validated within the workspace. Is this possible?


texts: Styling currently used feature of the Coverage layer of Atlas

texts:
I am looking for a way to define a different style just for the currently used feature of the Coverage layer of Atlas.
The image below will hopefully make it clearer: the left side is how my composition is looking, and on the right is what I would like to achieve. In my case the coverage layer is a map of land plots and I want to have the Page number feature to be styled differently to highlight it as the main object of the map composition.

I have managed to use data overrides to have styles that are applied only in the context of a Print Layout, but they all seem to operate at the layer level and not at a feature level.
Any ideas? Thanks


texts: arcpy.JSONToFeatures_conversion and json module conflict

texts:
I'm using arcpy.JSONToFeatures_conversion() to batch convert some json files to features in a gdb. Python kept running into errors (either ERROR 001558: Error parsing json file, or ERROR 999998). This was troubling since I had previously tested the input/output without issue. After a lot of troubleshooting I found the error only occurs after the json module was imported. What is the most pythonic way of dealing with this? Can anyone explain the conflict and how I might recognize/anticipate this when troubleshooting code in the future?

Full error:

Traceback (most recent call last):
  File "<pyshell#10>", line 1, in <module>
    arcpy.JSONToFeatures_conversion(in_json_file, out_features)
  File "C:\ArcGIS\Desktop10.3\ArcPy\arcpy\conversion.py", line 528, in JSONToFeatures
    raise e
ExecuteError: ERROR 001558: Error parsing json file 'L:/<dir>results16870.json'.
Failed to execute (JSONToFeatures).

Code Snippet (works when json isn't imported):

import arcpy
import json
temp_json = r"~\result1216.json"
fcName = r"~\Testing.gdb\newFC"
arcpy.JSONToFeatures_conversion(temp_json, fcName)

Screen shot of code/process and error (note the same input/output):


texts: Creating tool for finding catchments, watersheds and channel network in ArcGIS for Desktop?

texts:
I've been working in QGIS for a while, but I have to use ArcGIS from now on. In QGIS I made a tool in the Graphical Modeler, in order to ease a process I had to do several times. I've been told this is possible in ArcGIS as well, but I would like some hints, as the names aren't the same as I'm used to.

My old tool used a DEM/Raster and then did the following (Using SAGA)

Fill Sinks -> Catchment Area -> Channel Network -> Watershed Basins

Simple tool for finding catchments, watersheds and the channel network. Which tools would do something similar in ArcGIS? I've been reading up on it, but there seems to be several?

The function of the tool is finding where the water flows, and the catchments/watersheds to these channels.


texts: Adding input files of near tool using %name% or %value% function?

texts:
I am trying to caluclate the nearest distance between road and road_buf shp.

But I'm struggling to add input features using %name% function in NEAR TOOL.

I want to add road_buf file as "c:Desktop\road_BF\%name%_bf.shp",,,
but it does not work.

Does not NEAR TOOL provide %name% function as input ?


texts:
I have an Openlayers map set up, it is all working well except in IE7 (which most of our users have) where performance is slow.
I load in several layers of vector features and then switch these layers on/off using a custom 'layer switcher' which is some buttons/icons on the webpage.

In Chrome this works well, it's smooth and fast. IE7 is slow and there is usually a big delay when switching on a certain layer with 250 point features. The delay itself is ok, but it's easy for a user to think it hasn't worked and click the button again.

I decided to temporarily replace the button image with a loading image, make layer visible, then change button back to its correct image:

function toggle(buttonClicked) {
        button = document.getElementById(buttonClicked);
        var layer = map.getLayersByName(buttonClicked);
        if ( layer[0].visibility === false) {
            toggleButtonImg(buttonClicked, 'loading');
            layer[0].setVisibility(true);
            toggleButtonImg(buttonClicked, 'on');
        } else {
            layer[0].setVisibility(false);
            toggleButtonImg(buttonClicked, 'off');
        }

}

And the toggleButtonImg function:

function toggleButtonImg(buttonName, status){
        button = document.getElementById(buttonName);
        if (status === "loading") {
            button.src = "images/loading.gif";
        }else {
            button.src = "images/" + buttonName + "_" + status + ".png";
        }
    }

And if it helps, the button code:

<input type="image" onClick="toggle('Postcode')" id="Postcode" src="images/postcode_off.png">

The issue is, when you click the button the image doesn't change, the browser 'locks up' as it is loading the layer, then the button changes to the 'on' version. (off = greyed out, on = solid).
I've added some extra code to the toggleButtonImg() in the hope the image would update before loading the layer:

button.src = "images/custSymb/loading.gif";
alert("hey the button should now have changed");

but I click, image stays 'off', alert appears, I ok it, image is still 'off', layer loads and image changes to 'on'.

Any ideas on how I can get around this issue and also any other tips to improve IE performance.
e.g Will postcode_layer.removeAllFeatures(); be beneficial to perfomance when the postcode layer is not visible? (I realise I will have to reload it each time).

Thanks


texts: How should I go about calculating the centroid of several lat long points in python, and perl?

texts:
Calculating the centroid several time on a number of geographical points. I was wonder is there a perl or python module that could make my life easier when working with large datasets of lat long points?

A coworker recommended using the perl module Math::Polygon, and that the polygon_centroid() function could be of some use, however I would like to know of other options.

Thank you for your time,


texts: Naming convention for scheduled geoprocessing outputs

texts:
I'm planning on implementing several geoprocessing python scripts that are scheduled to run nightly. The scripts will either output new feature classes, or update existing ones. The environment is an Oracle geodatabase.

I find with other database objects, it's helpful to be able to easily differentiate between object types, sources etc. in ArcCatalog (or other database software). Examples would be:

Feature classes and non-spatial tables have different icons in the catalogue tree
Specialized geodatabase objects have unique icons
Views are suffixed with _VW (best practice)
Versioned views are suffixed with _EVW

However, when it comes to scheduled geoprocessing outputs, I'm not aware of any way to differentiate them from regular feature classes at first glance. Sure, metadata is an option, but I'm looking for something a little more in-your-face than that.

Is there any sort of naming convention that is used for scheduled outputs? Maybe suffixing the table with a _GP (for geoprocessing) or _NT (for nightly task)?

For example:

I have read-only access to a non-spatial INCIDENTS table from 3rd party software. The table has X and Y columns. I have a geoprocessing script (scheduled nightly task) that processes the table. It creates XY event layer, does a spatial join, and outputs minimum bounding geometry. In other words, the output is a copy of the original data; the only difference is that it has been processed.

I'm looking for a way to identify geoprocessed copies of data, such as the example mentioned above. It would help us administrators easily keep track of which tables are master tables, and which ones are copies.

It could be something as simple as a _COPY suffix. Or maybe there are established database best practices/conventions.


texts: In a country I have numeric data by province, but not for all the provinces. With arcgis I would like to interpolate those data to the no-data provinces. Of course I have administrative georeferenced shapefiles.


texts: How can I avoid UnicodeEncodeError when using geoprocessing tools in QGIS?

texts:
I have a shapefile, that I want to split with the "Split Layers by Field"-tool. Due to some umlauts (ä,ü,ö) there will be a python error as soon, as the tool reaches a feature with an umlaut. The following error message will appear:

  *An error has occured while executing Python code:
  Traceback (most recent call last):
    File "C:/Users/Gidi/.qgis//python/plugins\layers_by_field\layers_by_field_dialog.py", line 58, in accept
      self.split( inLayer, inField )
    File "C:/Users/Gidi/.qgis//python/plugins\layers_by_field\layers_by_field_dialog.py", line 146, in split
      self.vlayer = QgsVectorLayer(vProvider.dataSourceUri(), str(layer.name()) + "_" + str(uValues[j]), "ogr")
  UnicodeEncodeError: 'ascii' codec can't encode character u'\xfc' in position 7: ordinal not in range(128)
  Python version:
  2.7.2 (default, Jun 12 2011, 15:08:59) [MSC v.1500 32 bit (Intel)]
  QGIS version:
  1.8.0-Lisboa Lisboa, 6416f38
  Python path: ['C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python', 'C:/Users/Gidi/.qgis//python', 'C:/Users/Gidi/.qgis//python/plugins', 'C:/PROGRA~2/Quantum GIS Lisboa/apps/qgis/./python/plugins', 'C:\PROGRA~2\Quantum GIS Lisboa\bin\python27.zip', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\DLLs', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\plat-win', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\lib-tk', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\bin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\PIL', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\win32\lib', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\Pythonwin', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\Python27\lib\site-packages\wx-2.8-msw-unicode', 'C:\PROGRA~2\Quantum GIS Lisboa\apps\qgis\python\plugins\fTools\tools']*

What can I do to avoid this error? I do not want to change the name of every feature and would prefer the umlauts to remain in my shapefile.


texts: How to count number of stores per company per region?

texts:
I have two layers. One point layer and a Polygon layer.
It is simple to count the points in the polygon layer by the "points in polygon" function. But I also need some more information of the point layer in the polygon layer.

For example:

The point layer includes information of the company name (and some more) of every point.

--> I need the number of Points in the polygon layer but also some information about how many stores of company XYZ are in the polygon

I tried my luck with the "select by location" function but this works just the other way round (I can create a point layer with all necessary information but not a polygon layer because the function only uses the "attributes of the first exact object".

The other possibility of chosing a summarization is not in my interest because I need to keep all information - anyway it does not work because my column names will be to long after creating it.

I will add a screenshot to better explain my problems (it is the german version and I am actually working in France so please try to overlook my faults ;)


texts: What is the structure of Standard Interchange Format (SIF) used by spatial data providers?

texts: What is the structure of Standard Interchange Format(SIF) used by navteq and other spatial data providers? Any documentation/sample files? Is there a SIF parser available in C/C++/Java/PHP ?


texts: Using just the variable name (not full file path) for inline variable substitution in Model Builder

texts:
I am trying to use inline variable substitution in ModelBuilder to return the name of a raster used in Zonal Statistics as a field in my output table (see model below).

Currently, the model works fine using the "%Raster Dataset%" to get the field I want, but it returns the FULL file path - e.g. C/Users/Newton/GIS/Project/RASTER123 instead of just "RASTER123". Anyone know how to just select the name of the input raster, not the full path?


texts: Eliminating/minimizing distortion in 200 NM EEZ buffers?

texts:
I frequently need to calculate buffers of 200 nautical miles from country or baseline polygons in order to map maritime Exclusive Economic Zones (EEZs). However, I've noticed that when generating the buffers in QGIS, I get different results depending on the projection specs of the layer I'm buffering from. For example, in a recent project I'm using an azimuthal equidistant projection, and the buffer turns out differently depending on where the projection is centered. This also means that even if I center the projection on the country, there is often still significant distortion at more distant corners of the same country.

Is there any way to generate buffers without any distortion (i.e. equidistant from each node on the polygon from which they're being calculated)?


texts: How to conditionally format a table from Excel in the layout view?

texts:
I'm making a layout view, thing and have run into a problem. I have formatted my table beautifully in Excel, see Picture 1, but when I import it to QGIS it clears the formatting. The color goes away, and more annoyingly my percentages become 10 point decimals, oh no, Picture 2!
How can I implement the table like you see in the first in the print layout element. I know I could just insert the picture of the Excel document, but I want to know how to do this in theory with a attribute table made in QGIS.


texts:
I am trying to display certain features in a shapefile-based layer using QGIS. I have found the subset function in the properties panel, which does exactly what i want, but is too tedious given that I have to manually sift through a large number of features. Is there a quicker way to do this?

I am using both shapefiles and postGIS
best,
dietmar


texts: How to extend Approximal Medial Axis both sides

texts:
I would like to extend lines created with the function ST_ApproximalMedial_Axis to 2 points, located each side of this line.
Each line is the medial axis of a polygon created by the intersection of 2 buffers and the points are the points of intersection of these buffers.

    -- Creation of axis
    CREATE TABLE public.l_axe_med AS
    SELECT
    row_number() OVER () as id_axe,
    t3.idu1,
    t3.idu2,
    ST_ApproximateMedialAxis(t3.geom) geom
    FROM (
    SELECT t1.idu as idu1, t2.idu as idu2, ST_INTERSECTION(t1.geom, t2.geom) geom
    FROM public.tb_buffer_initial t1
    INNER JOIN public.tb_buffer_initial t2
    ON ST_INTERSECTS(t1.geom, t2.geom) IS TRUE AND t1.geom && t2.geom
    WHERE t1.id_ilot != t2.id_ilot and t1.idu1 != t2.idu2) t3

    -- Creation of the points of intersection

    CREATE TABLE public.p_interbuf AS SELECT
    row_number() OVER () as id_interbuf,
    t1.id_buf as id_buf1,
    t2.id_buf as id_buf2,
    t1.idu as idu1,
    t2.idu as idu2,
    ST_INTERSECTION(ST_BOUNDARY(t1.geom),ST_BOUNDARY(t2.geom)) geom,
    ST_GeometryType(ST_INTERSECTION(ST_BOUNDARY(t1.geom),ST_BOUNDARY(t2.geom))) as type_geom
    FROM
    public.tb_buffer_v2 t1
    INNER JOIN public.tb_buffer_v2 t2
    ON ST_INTERSECTS(t1.geom, t2.geom) IS TRUE AND t1.geom && t2.geom
    WHERE t1.id_ilot != t2.id_ilot and t1.idu != t2.idu and t1.id_buf != t2.id_buf

I imagine there will be a call of the function ST_Closest... but I have for the moment no idea of how to proceed. May you help me please ?

I've found a part of the answer with the following request

SELECT
id_axe,
idu1,
idu2,
CASE WHEN ST_Distance(ST_StartPoint(axe_geom), geomp1) < ST_Distance(ST_StartPoint(axe_geom), geomp2)
THEN ST_MakeLine(array[ST_MakeLine(geomp1, ST_StartPoint(axe_geom)), axe_geom, ST_MakeLine(ST_EndPoint(axe_geom), geomp2)])
ELSE ST_MakeLine(array[geomp2, axe_geom, geomp1]) END as geom
FROM (
SELECT
id_axe,
idu1,
idu2,
axe_geom,
ST_GeometryN(multipoint_geom, 1) as geomp1,
ST_GeometryN(multipoint_geom, 2) as geomp2,
ST_NPoints(t4.multipoint_geom) as nb_path
FROM
(SELECT
row_number() OVER () as id_axe,
t3.idu1,
t3.idu2,
ST_ApproximateMedialAxis(t3.poly_geom) as axe_geom,
t3.multipoint_geom
FROM (
SELECT t1.idu as idu1, t2.idu as idu2, ST_INTERSECTION(t1.geom, t2.geom) as poly_geom, ST_INTERSECTION(ST_BOUNDARY(t1.geom), ST_BOUNDARY(t2.geom)) as multipoint_geom
FROM public.tb_buffer_initial t1
INNER JOIN public.tb_buffer_initial t2
ON ST_INTERSECTS(t1.geom, t2.geom) IS TRUE AND t1.geom && t2.geom
WHERE t1.id_ilot != t2.id_ilot and t1.idu != t2.idu) t3) t4) t5

But whatever the case, the line that is created is a the line between the intersections of boundaries although I Would like to use the true geometry of my medial axis.

And when I use ST_Union I only get back the medial axis, but it is not extended to the boundaries

THEN ST_Union(array[ST_MakeLine(geomp1, ST_StartPoint(axe_geom)), axe_geom, ST_MakeLine(ST_EndPoint(axe_geom), geomp2)])
ELSE ST_Union(array[ST_MakeLine(geomp2, ST_StartPoint(axe_geom)), axe_geom, ST_MakeLine(ST_EndPoint(axe_geom), geomp1)]) END as geom


texts: Get lat/long given current point, distance and bearing

texts:
While searching the web, I've encountered the following formulas for calculating the destination point P2's location given the source point P1 = (lat1,long1), the distance from it & the bearing:

lat2 = math.asin( math.sin(lat1)*math.cos(d/R) +
   math.cos(lat1)*math.sin(d/R)*math.cos(brng))

lon2 = lon1 + math.atan2(math.sin(brng)*math.sin(d/R)*math.cos(lat1),
   math.cos(d/R)-math.sin(lat1)*math.sin(lat2))

Where P2=(lat2,lon2).

I was trying to develop the math behind it but could not really get it right

Can someone please show me how to develop the lat2 & lon2 mathematically?


texts: Using TIN - interpolation plugin - to calculate volume

texts:
I am having trouble calculating the volume of my TIN.

I have attempted to create a TIN using the interpolation plugin in QGIS. Here are the steps that I followed.

I added a Zvariable row in the attribute table of the vector layer with the value of 500 (I am assuming this corresponds to 500 meters).
I used the interpolation plugin in QGIS to create a TIN layer. I created it following steps advised through the qgis website which recommended however that I use a cell size of 5000 for both both X and Y cells but this came up with an error message for me so I left it as it was set by default (please see the screenshot below for details of exactly what I entered).

After, I defined the coordinate system as WGS84, and the attached layer appeared (please see in the below attachment, the layer titled testTIN which is shown in black shading).

Do you believe that my TIN layer is correct and if so, how do I now go about calculating the volume of the area?

I have also attached for you a screenshot of the attribute table of the layer used in the interpolation plugin to create the TIN.


texts: .MXD file shows as *.DOC file after upgrading machine with ArcGIS 10.3.1 for Desktop to Windows 10?

texts:
I created an .mxd file on a computer with an earlier version of Windows (I was at a college campus; not sure what version). Now, on a computer with Windows 10, the file appears as a Word .DOC and ArcMap cannot open it.  My only idea was to reassign the default software type that opens the file, but that did not work.

My mxd is a word file? cites a similar problem and says installing ArcGIS 10.3.1 on the Windows 10 machine will fix it, but I already have ArcGIS 10.3.1 installed. Also I don't think that would be the problem, because I can open other .mxd files, and Windows Explorer is interpreting the file as .DOC regardless of ArcGIS being installed or not.

Any tips for converting this back to an .mxd I can open on a Windows 10 machine using ArcGIS 10.3.1?


texts: Color region polygons according to cluster label?

texts:
I want to present data on a map of our country.

I have used R for showing the clusters on a country! Here the steps I took:

setwd('D:/r/cluster2')
channel <- odbcConnectExcel('cluster.xls')
data <- sqlFetch(channel, 'clust9')
y9 <- data.frame(inf=data$infest, faible=data$faible, moyen=data$moyen, fort=data$fort, lon=data$Lon, lat=data$Lat)

library(fossil)

d = earth.dist(y9)
km <- kmeans(d,centers=5)
hc <- hclust(d)
clust <- cutree(hc, k=5 )
set.seed(123)
plot(hc)

y9$clust <- cutree(hc,k=5)
map.AL  <- readOGR(dsn="D:/r/cluster/shp", layer="ALG_boundaries")
map.df  <- fortify(map.AL)

ggplot(map.df)+ geom_path(aes(x=long, y=lat, group=group))+
  geom_point(data=y9, aes(x=long, y=lat, color=factor(clust)), size=4)+
  scale_color_discrete("Cluster") + coord_fixed()

On the map, it's showing just points for the regions that have the same characteristics. However, what I want to show is all the region colored not just a point.


texts: Calculate field with a value in a list based on field

texts:
I don't want to write an if statement for every field. I want to know if there is another way to field calculate each field with the values related to it in this example table:

Fields
(X,Y,INSTALLDATE,BLDGID,LIFECYCLE,FACILITYID,REGION)

Values
(0,0,1/1/9999,NA,TBD,NA,NA)

import arcpy
arcpy.env.overwriteOutput = True
arcpy.env.workspace = Database Path
Asset = "Asset"
try:
    for fc in arcpy.ListFeatureClasses("ONUS_ssManhole") + arcpy.ListFeatureClasses("ONUS_ssInterceptor") + arcpy.ListFeatureClasses("ONUS_ssHoldingTank") + arcpy.ListFeatureClasses("ONUS_ssCleanOut") + arcpy.ListFeatureClasses("ONUS_ssFitting") + arcpy.ListFeatureClasses("ONUS_ssLiftStation") + arcpy.ListFeatureClasses("ONUS_ssValve") + arcpy.ListFeatureClasses("ONUS_ssPressurizedMain"):
        fcs = fc
        arcpy.MakeFeatureLayer_management(fc, Asset)
        for field in arcpy.ListFields(fcs, "X") + arcpy.ListFields(fcs, "Y") +  arcpy.ListFields(fcs, "INSTALLDATE") +  arcpy.ListFields(fcs, "LIFECYCLE") + arcpy.ListFields(fcs, "PROJECTID_ADD") +   arcpy.ListFields(fcs, "FACILITYID"):
            fields = field
            expnull = ("{0} is NULL").format(field.name)
            arcpy.SelectLayerByAttribute_management(Asset, "NEW_SELECTION", expnull)
            result = arcpy.GetCount_management(Asset)
            count = int(result.getOutput(0))
            if count == 0:
                print(fc + " {0} is Good").format(field.name)
            if count > 0:
                print (count)
                arcpy.CalculateField_management(Asset,(????????)
                print(fc + " {0} is now clear of NULL Values").format(field.name)
except Exception:
    e = sys.exc_info()[1]
    print(e.args[0])


texts: Need pros and cons of performing a GeoTIFF from ERDAS img conversion

texts:
I have a 20GB ERDAS IMG file that I need to incorporate into my project.
It's data covers the entire USA, indicating the thickness of satellite sensed tree canopy.
What image format would provide me with far better compressibility, (20 GB is a lot of data!) but still allow me to access differing tiled sections of the data without trying to process or load the whole image.
Or would it be better to just section the file off myself and attach it in tiled pieces to a database myself?
The tiling aspect of ERDAS Imagine makes it sound attractive for searching and indexing, but the size is alarming.


texts: Extract custom shape area from Openstreetmap pbf

texts:
I have a full-country data download of OSM from http://download.geofabrik.de/ . I want to extract OSM data for a target area that is thin and winding/irregular in shape. I could take the max-min lat-long bounds,

bounds="78.55834,12.54397,80.3558,13.71654"
osmconvert india-latest.osm.pbf -b=${bounds} --complete-ways -o=area1.pbf

.. But that leads to lot of wastage as the rectangle pbf created is many times heavier than the actual area needed. Actually, by taking full rectangle, some settlements with very dense OSM mapping are being included which balloon up the output pbf size like anything. My target area is mostly rural with sparse data. This data has to be churned by a process where minimizing pbf size is of essence to achieve best performance.

If I have my area as a shape in a standard .geojson format shapefile, how do I "clip" the pbf precisely by this shape?


texts:
I need the twisted framework (Python) for a qgis plugin on Windows 7.
I installed it with the pip installer.
When I want to start the plugin it doesn't find twisted.
When I start a python interpreter in the cmd.exe twisted can be imported.

this is my path variable: C:\Python27\Lib\site-packages;C:\BIN;C:\Python27\Scripts;C:\Python27;

twisted is installed inside C:\Python27\Lib\site-packages

How should I install a python package when I want to use it inside QGIS?


texts: R does not load coordinate reference when importing a multiband tif raster from Triplesat satellite

texts:
I'm trying to read Triplesat imagery into R. However there seems to be a problem with the coordinate reference in the .tif file as the crs, extent and resolution are not loaded. I tried two ways of loading the data with raster() from the raster package and readGDAL() from rgdal. Both have the same problem, other .tif files work fine. QGIS however does open the .tif and displays it properly because it finds the extents in a text file that is located in the same folder.

R information on raster object:

class       : RasterLayer
band        : 1  (of  4  bands)
dimensions  : 7775, 7810, 60722750  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 7810, 0, 7775  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : D:\\Triplesat\\MyRaster.tif
names       : TRIPLESAT_3_MS_L1_20170802093353_001097VI_006_20170630001001_030
values      : 0, 65535  (min, max)

My question is, how can I open the .tif in R so that the coordinate reference is properly loaded? Either by getting the values from a metadata-file or by defining the coordinate reference after that. If possible I would like to stay within R.

What I have tried:

Assigning a crs and extent with

r <- raster("D:\\Triplesat\\MyRaster.tif")
crs(r) <- "+init=epsg:4326"
extent(r) <- extent(reference_raster) # reference raster has the same spatial boundaries

This however does not work because the extent is based on the min.x, min.y, max.x, max.y values which creates a rectangle bounding box with the x-axis in east-west direction. The raster however is not so straight, see:

Loading the data in QGIS and exporting it so R can open it. This works for a few files but I need to perform the same procedure dozens of times so this method is too time consuming.
Using ESA's SNAP Desktop. However also SNAP does not load the coordinate references.

WeTransfer download link to files: https://we.tl/WdwNMhvqFH

Edit:

Included in the files is also a browser file that is used for quick viewing the data. This does load properly in R and gives me the correct crs information. Also in QGIS the projection is displayed in lattirude longitude. See the R result:

browser <- raster("D:\\Triplesat\\MyBrowserRaster.tif")
crs(browser)
CRS arguments:
     +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0


texts: Selecting WMS backdrop in QGIS print composer?

texts:
I have a WMS service from Geoserver which is controlled by a parameter (SQL params). This means that the following string is added to the getmap-request to Geoserver:

viewparams=objectid:11

With this I will get a tailored wms-map for a category of objects. I am thinking of creating an atlas based on the a control table which holds information about the objectid, plate name, description and more.

platenumber  objectid  objectname   description
1            10        Polar bears  This map shows where to find polar bears.
2            20        Crows        This map shows where to find crows.
3            55        Dragons      This map shows where to find dragons.

Is it possible to dynamically (through the atlas) control the wms getmap string?


texts: Is it possible to underline text in a Google Earth Engine panel?

texts:
I am currently writing a code for an app in which I would like to underline some text. However, looking at the ui.Panel.style options there is none which allows underlining the text.
It says that style properties behave like their CSS counterparts, but the CSS property for underlining is
text-decoration: underline

But this returns an error on GEE (Invalid parameter for layout manager). I tried textDecoration just in case but it didn't work either. Is it even possible to underline text? Or is this a dead end?
Here is an example:
 var panel = ui.Panel();
  panel.style().set({
    width: '25%',
    position: 'bottom-right',
    maxHeight: '80%'
  });

Map.add(panel);

var linkLabel = ui.Label('Example',{fontSize: '12.5px'/* other properties go here*/}).setUrl('http://www.example.com/');
var linkPanel = ui.Panel([linkLabel]);

panel.add(linkPanel);


texts: Cannot get L.control.coordinate to work properly

texts:
I am trying to make a map modal with Leaflets that displays the lat./lng. coordinates in a box when the user clicks on the desired point. I got the popup given in the Leaflets tutorial to work fine and now I am trying to use the Controls.Coordinates plugin. I currently using the the Coordinates Control plugin, but I keep getting errors like NumberFormatter.js:1 Uncaught ReferenceError: L is not defined (Note: I have the Control.Coordinates.js installed) or Uncaught TypeError: L.Control.Coordinates is not a constructor. I think there must be a problem with how I have my plugin path or perhaps the javascript. Any ideas?

<title>Create Job</title>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="stylesheet"
href="http://maxcdn.bootstrapcdn.com/bootstrap/3.3.6/css/bootstrap.min.css">
<link rel="stylesheet" type="text/css" href="./main.css">
<link rel="stylesheet"
href="http://cdn.leafletjs.com/leaflet/v0.7.7/leaflet.css" />
<link href='https://api.mapbox.com/mapbox-gl-js/v0.19.1/mapbox-gl.css'
rel='stylesheet' />

<link rel="stylesheet" href="./src/Control.Coordinates.css"/>

And here is the javascript part. The problem seems to occur when I declare var c = L.Control.Coordinates():


texts: MapServer WMS background color in OpenLayers 3

texts:
I am using MapServer to render a number of layers and then I'm setting it up as a WMS server in order to render it using OpenLayers 3. I want the background color of the image to be grey so I set it with the following snippet:

MAP
 EXTENT 300000.000 4374000.000 316000.000 4386000.000
  UNITS meters
  SIZE 900 580
  IMAGETYPE png
  IMAGECOLOR 240 237 229
  ...

So in CGI mode it behaves as expected rendering the background as grey. However the same does not happen with OpenLayers. I'm setting params to render it as 'image/png':

var params = {
            LAYERS : 'layer1, layer2,layer3...layerN',
            FORMAT : 'image/png'
     }

Setting FORMAT as 'image/png' makes the background color white instead of the grey I want. If I set FORMAT as 'image/jpeg' the the background color turns to total black. So how can I keep the same grey color using OpenLayers?

Which format do I need to use?


texts:
I am trying to execute a buffer operation using geoserver 2.8.3, but i'd rather use json than GML to send this request. Is it possible to perform an Execute operation this way?

The describe process request contains this little snippet for data input

            <Input maxOccurs="1" minOccurs="1">
            <ows:Identifier>geom</ows:Identifier>
            <ows:Title>geom</ows:Title>
            <ows:Abstract>Input geometry</ows:Abstract>
            <ComplexData>
                <Default>
                    <Format>
                        <MimeType>text/xml; subtype=gml/3.1.1</MimeType>
                    </Format>
                </Default>
                <Supported>
                    <Format>
                        <MimeType>text/xml; subtype=gml/3.1.1</MimeType>
                    </Format>
                    <Format>
                        <MimeType>text/xml; subtype=gml/2.1.2</MimeType>
                    </Format>
                    <Format>
                        <MimeType>application/wkt</MimeType>
                    </Format>
                    <Format>
                        <MimeType>application/json</MimeType>
                    </Format>
                    <Format>
                        <MimeType>application/gml-3.1.1</MimeType>
                    </Format>
                    <Format>
                        <MimeType>application/gml-2.1.2</MimeType>
                    </Format>
                </Supported>
            </ComplexData>
        </Input>

which gives me some hope that a full json request might be possible but it currently doesn't work when i try it, and i can't seem to find any working examples out there.

The json data i used unsuccessfully is

{
    "wps:Execute": {
    "-version": "1.0.0",
        "-service": "WPS",
        "-xmlns:xsi": "http://www.w3.org/2001/XMLSchema-instance",
        "-xmlns": "http://www.opengis.net/wps/1.0.0",
        "-xmlns:wfs": "http://www.opengis.net/wfs",
        "-xmlns:wps": "http://www.opengis.net/wps/1.0.0",
        "-xmlns:ows": "http://www.opengis.net/ows/1.1",
        "-xmlns:gml": "http://www.opengis.net/gml",
        "-xmlns:ogc": "http://www.opengis.net/ogc",
        "-xmlns:wcs": "http://www.opengis.net/wcs/1.1.1",
        "-xmlns:xlink": "http://www.w3.org/1999/xlink",
        "-xsi:schemaLocation": "http://www.opengis.net/wps/1.0.0 ../../xsd/wpsAll.xsd",
        "ows:Identifier": "JTS:buffer",
        "wps:DataInputs": {
        "wps:Input": [
            {
                "ows:Identifier": "geom",
                "wps:Data": {
                    "wps:ComplexData": {
                        "-mimeType": "text/xml; subtype=gml/3.1.1",
                        "#cdata-section": "POINT(-9128198.245048434 5311971.846945471)"
                    }
                }
            },
            {
                "ows:Identifier": "distance",
                "wps:Data": { "wps:LiteralData": "1000" }
            },
            {
                "ows:Identifier": "quadrantSegments",
                "wps:Data": { "wps:LiteralData": "0" }
            },
            {
                "ows:Identifier": "capStyle",
                "wps:Data": { "wps:LiteralData": "flat" }
            }
        ]
    },
    "wps:ResponseForm": {
        "wps:RawDataOutput": {
            "-mimeType": "text/xml; subtype=gml/3.1.1",
                "ows:Identifier": "result"
            }
        }
    }
}


texts:
I'm using Carto Mobile in an Android app and I wonder how could I align the Carto logo on a corner that I choose, instead of having it randomly displayed and interfere with map controls I'm displaying over the map.

Usually it appears in the bottom left corner, and I have a scale bar over there. When it gets displayed there, I cannot see the scale bar correctly. I'd prefer the Carto logo to always be in the bottom right corner, since I have nothing over there.

Aside from that, I guess that if I get to pay for Carto, the logo could be removed, isn't it?


texts:
I am trying to work out how much bigger is Coverage1 over Coverage2 in Mapinfo pro.

I have 2 Coverage data, Coverage1 and Coverage2.  Both maps are in .tab file format.  Coverage1 has a wider footprint than Coverage2.  Coverage1 and Coverage2 also overlaps in places.  I would also like to point out that, Coverage data is scattered.  I am looking to perform Spatial query to highlight the additional coverage provided by Coverage1, essentially (Coverage1 – Coverage2).

I am currently exploring the “Combine functionality” in Mapinfo pro as a solution.


texts: Increase speed of reverse clip (erase) for large shapefiles in R

texts:
I would like to do a reverse clip (or erase in ArcMap) from shapefiles (polygons) then to merge the two shapefiles into a single shapefile. From this post  How to 'reverse clip' (erase) in R, I tested the function gDifference from the rgeos package but it's very long (several days). How can I improve the speed or is there another way to do a reverse clip in R?

Here are some information on the two shapefiles (I kept the ArcMap terminology http://pro.arcgis.com/en/pro-app/tool-reference/analysis/erase.htm)

## Input feature
> shp1
class       : SpatialPolygonsDataFrame
features    : 61378
extent      : -74, -71.5, 45, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
variables   : 55

## Erase feature
> shp2
class       : SpatialPolygonsDataFrame
features    : 429418
extent      : -79.53649, -61.75487, 44.99248, 50.25604  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0
variables   : 4


texts:
I have a MapDocument that I am trying to export to PDF. I am using this example. However, it always says it is finished but when I open the pdf file there is nothing in it. I have tried exporting to a TIFF but with the same results. Is there a known reason why this might be happening?

When it brings up the ArcMap program, it displays everything properly and I can manually export to PDF from the program just fine.


texts: Unable to get gdal_polygonize function to work in R

texts:
I know there have been previous questions with respect to the same. I have been trying to polygonize a raster using modified code from John Baumgartner.

Link : https://gist.github.com/johnbaums/26e8091f082f2b3dd279

I have installed Python 2.7 as well as installed the Anacondas framework for writing code. In addition, I installed the latest version of OSGeo4w for Windows, but I keep getting this error.

> p <- polygonizer(r)
  Error in ogrInfo(dsn = dsn, layer = layer, encoding = encoding, use_iconv
   = use_iconv,  :
     Cannot open data source
      In addition: Warning message:
     running command '"python" "C:\Users\vr235\DOWNLO~1\gdal222\GDAL-2~1.2\swig\python\scripts\GDAL_P~2.PY"
     "C:\Users\vr235\AppData\Local\Temp\RtmpeGRpnN\file33e8d1c7717.tif" -q  -f "ESRI Shapefile"
      "C:\Users\vr235\AppData\Local\Temp\RtmpeGRpnN\file33e825205ec7.shp"' had status 127

Here, r is the raster I would like to polygonize. I am using this code as rastertoPolygons takes forever to finish running.

Any suggestions?


texts: MODIS HDF to Corrected Reflectance Bands 7-2-1. Change Reflectance

texts:
Convert HDF to Tiff easily. But I can not get the images to 7-2-1 bands.

Like this :

How can I get this image from HDF and what product should I use (MYD09GQ and etc)? Maybe using ArcGIS?


texts:
I have an ungeoreferenced .tif of which I know the spatial information. However, when I run the following command it looses all unique values. All pixels have the value 32, which is the upperleft value.

gdalwarp -t_srs EPSG:4326 -to SRC_METHOD=NO_GEOTRANSFORM -te 5.4990372 46.2904918 15.5783806 55.7887200 -tr 0.020038456088761 0.020038456088761 -r near -overwrite %s %s" % (inputfiles, outputfile)

What am I missing?


texts:
I have 70 raster landsat images with srid: 32647,32648,32649 and want to clip raster to polygon (shapefile, srid: 4326) in postgis.

I imported all raster into 1 table with parameter -s : 32648, transform srid of shapefile to 32648 but the result is not as I expected ( only overlap 40% the polygon)

So How can I fix that problem ?


texts:
I'm wondering if it is possible to create stacked line profiles, to display either geophysical data collected along a line, or geochemical data also along a line?

Datamine Discover has a feature that allows you to use point data, in this case either geophysical or geochemical data, create a line profile and then display it on the map.

The advantage is it can highlight anomalies as the data lines cross a feature at ~90 degrees.

The Datamine tool allows you to create multiple lines based on multiple attributes, scale the lines, and highlight lines over certain cutoffs (as seen with the blue 'fill' in the below figure).

Obviously Datamine's solution is fairly polished as it is their business, but I was wondering if there was a workaround to have a similar coloured line effect in QGIS?


texts: How many forest clear cuts and how big are they?

texts: I have forest stand characteristic data from northern Sweden. We have interpreted multiple years of aerial photos and now have polygons outlining different forest stand ages over time (e.g. clear cuts, young forests, old forest, etc.). I'm interested in determining how the number of clear cuts has changed over time and how much area they have covered over time. But because of the way we mapped the forest stands, we can have 2 or more polygons that are called "clear cut" right next to each other. Thus, if I were to just count the number of polygons that were coded as "clear cut" this would over estimate the number of clear cuts. How can I get GIS (I'm a novice, but use both ArcGIS and QGIS) to count polygons that neighbor each other as just one polygon?


texts: Seeking method for determining survey points?

texts:
I'm trying to write software to create a topographic map, my goal is to select several points of interest and create a minimal mapping that is inclusive of the points and interpolates the topography.

Is there a standard algorithm / approach for progressively selecting additional points?

Update (1): @JoshC

My goal is to produce a mesh, but I'm not as interested in the x,y offset being consistent.

the simple case is with 2 points of interest

Let's say we have a point at (2,2,0) [elevation of 0] and another point of interest at (4,4,8).

The parameters of my algorithm will determine that a z delta of 8 units over a distance of 5.6 units is 'significant' and warrants additional exploration. I'll use the midpoint for this (4,4,?). So the software will recommend that I find z @ (4,4). Let's say I find that to be (4,4,2). From here, my algorithm will decide that a z delta of 2 over the distance of 2.8 [(2,2,0),(4,4,2)] is 'insignificant' but that the 6 unit increase over a distance of 2.8 [(4,4,2),(6,6,8)] is 'significant' and warrants additional exploration ie. surveying (5,5,?)

I'm satisfied with assuming radial topography, by minimal mapping I mean that I'll 'zoom' to a window of [(2,2),(6,6)]

The black box is the viewing window.

My thought is to then 'fill' the colored spaces with the weighted average of the z value and it's distance from the corresponding survey markers.

What I've described above is the algorithm I'm thinking about implementing, but I'm curious if there is already an industry standard (this is generalized, the solution needs to allow for an arbitrary number of points of interest / not just 2)

I'm new to this domain.


texts: Buffer Polygons change shape after running Dissolve Tool

texts: I am trying to dissolve my buffer polygons but when I do so the shape changes while dissolving. See pictures below. I am starting with point features, then I buffer them afterwards I want to dissolve the buffer polygons. I am using ArcMap 10.3. Any ideas on how to fix this issue?


texts:
I drew a line between two points and would like to change the line styles but currently is not working as expected. Tried to add the style into the features. Am i missing something?

        var points = [
            new OpenLayers.Geometry.Point( lng, lat ).transform(epsg4326, projectTo),
            new OpenLayers.Geometry.Point( qlng, qlat  ).transform(epsg4326, projectTo)
        ];

        var styles = {
          strokeColor: '#000000',
          strokeOpacity: 2,
          strokeWidth: 10
        };

        feature = new OpenLayers.Feature.Vector(
                new OpenLayers.Geometry.LineString(points, null, styles)
        );

        vectorLayer.addFeatures(feature);


texts: How to UPDATE a column with cluster IDs in PostGIS?

texts:
I have a table with points in it, and I want to cluster them.  I had created a separate table to hold the clusters, but I ended up copying all the data from the original table into the new cluster table.  Instead, I would like to add a cluster id (cid) column to the original table, and update it with the cluster id.  I tried

INSERT INTO spatial_data (cid)
SELECT st_clusterdbscan(geog::geometry, eps := 0.01, minPoints := 5) over() as cid
FROM spatial_data;

but it added a lot of new rows that were empty except for the cluster id.

I suspect I need to use the update command instead, but I don't know how to write the WHERE clause.

What is the correct way to populate the new cluster id column?


texts: Why is my Geoserver map off to one side?

texts: I am trying to load a map to display search results and when the map is loaded, it does not fill the modal that is supposed to contain it. I have tried changing the CRS between epsg:4326 and epsg:3857 to see if it was the projection that was messed up. Here is a screenshot of the problem, it should fill the entire modal:


texts: My client requires weather data for an application that has wind speed, wind direction aggregated by zip code for individual days of the week.  Where can this be acquired to be used in an Esri Web Mapping Application?  I have the capabilities to integrate this for a custom widget solution.


texts: Sum the values of all rasters in a user defined list

texts:
I am creating a QGIS plugin a part of which is like the SAGA Rasters sum tool. The user interface will have a widget to which the user will load the raster layers. I am trying to get the pixel-wise sum of values of the loaded rasters.

Like if the user loads 3 rasters R1, R2, R3 then for a pixel Pij, having values 5, 7 and 6 in R1, R2 and R3 respectively, in the new raster Pij = R1ij+R2ij+R3ij = 5+7+6 = 18

I am stuck up with both - creating the UI and the calculations.

So the first part of my question is about creating a QListWIdget and dragging and dropping the layers into the widget and using these layers for additions. The best I could get to was to add the QListWidget into the UI using QtCreator. Couldn't get to add and call layers from the QListWidget.

Here is the code I used with a push button to add layers in the QListWidget, but the layers add as a single row and I don't know how to use these for the actual calculations.

def add_rasters(self):
    name=QFileDialog.getOpenFileNames(self.dlg, 'Select Rasters', '','*.tif')
    name_str = str(name)
    list = []
    list.append(name)
    self.dlg.listWidget.insertItem(2,name_str)

self.dlg.pushbutton.clicked.connect(self.add_rasters)

Using the above code the ListWidget gets filled up like:
(['C:/QGIS_plugin/rstr1.tif', 'C:/QGIS_plugin/rstr2.tif', 'C:/QGIS_plugin/rstr3.tif'], '*.tif')

Next I am querying about the calculations:

I can do the addition calculation by manually adding up all the layers in the list; but I could not figure out a way to do the calculation automatically for n-number of the layers loaded in the widget.

Below is my code for addition of layers already loaded in the QGIS map project.

import processing
from qgis.core import *
from PyQt5.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry

entries = []

myList=[layer for layer in QgsProject.instance().mapLayers().values()]
print(my_list)

for i in range(0, len(myList)):
     raster = myList[i]
     #readRst = QgsProcessingUtils.mapLayerFromString(raster)
     ras1 = QgsRasterCalculatorEntry()
     ras1.raster = raster
     ras1.ref = "ext_lyr" + str(i+1) +"@1"
     ras1.bandNumber = 1
     entries.append(ras1)

formula = '(' + entries[0].ref + ' + ' + entries[1].ref +  ')'

output = r"C:/QGIS_plugin/rstrs_sum.tif"
calc = QgsRasterCalculator(formula, output, 'GTiff', raster.extent(), raster.width(), raster.height(), entries)
calc.processCalculation()

I am afraid that I am unable to frame this question properly (mixing two questions), but I am myself not sure how this would work best. I have gone through several posts on this site and now I am completely lost!

I am quite new to python programming and QGIS, hence need help with the simplest things also. Hoping to get some help here with this.


texts: Exporting Origin-Destination Matrix Table into Separate Origin and Destination Columns

texts:
Trying to work with exported data tables from Network Analyst of an origin-destination matrix for further network analysis like centrality measures...etc., but having to re-organize the origin and destination into their separate columns for other software like Python or NodeXL to perform well.

Is taking up quite a bit of time to do manually, does anyone know of any workarounds?


texts: ST_Voronoi is creating invalid geometries

texts:
I have a workflow where I am creating Voronoi polygons from a series of points. However since the polygons on the exterior extend really far out I have clipped them and then joined them back to points to get the IDs back. However when I join them back based on a ST_Intersects about 1-2% of them get lost. Here is the original Voronoi polygons and then after the intersection.

Here is the code that I am using

CREATE OR REPLACE VIEW voronoi_test AS
WITH voronoi_temp2 AS (
    SELECT (ST_dump(ST_VoronoiPolygons(ST_Collect(the_geom)))).geom as the_geom
    FROM ways_vertices_pgr)
SELECT id, voronoi_temp2.the_geom
FROM voronoi_temp2
JOIN ways_vertices_pgr
ON ST_Intersects(ways_vertices_pgr.the_geom, voronoi_temp2.the_geom);

Upon further inspection of the geometries those that drop out have duplicated points. So something like [(0, 0), (1, 0.5), (1, 0.5), (0, 1), (0, 0)]. So they are self-intersecting and invalid. I have tried MakeValid and RemoveRepeatedPoints but both end up dropping the invalid polygons.

This only occurs when the data is in a geographic coordinate system and not when in a projected coordinate system.


texts: Creating button on Python AddIn toolbar to run python script when clicked?

texts:
I need to create a button on toolbar of arcgis 10.2 who run a python 2.7 script created before when I click. Then I need to know if it's possible in this button that script write results to an excel file.
To create a button I went to: http://desktop.arcgis.com/en/arcmap/10.3/guide-books/python-addins/creating-an-add-in-button.htm an I had follow steps, but I didn't know if this resulted... In these steps I never mentioned the document or the folder where I want this button will appear..


texts: Any QGIS plugin to facilitate the creation of hiking maps?

texts:
There used to be this one https://plugins.qgis.org/plugins/AreaPrinter-master/

whose description reads::

  AreaPrinter simplifies the process of creating print composers for
  navigation-quality hiking maps. It produces a print composer with km
  grid, scale and ruler for an area that may span several pages.
  Currently A4 paper, 1:10k,1:25k or 1:50k scales, and WGS 84 / UTM
  projections are supported. Warning: Verify the correctness of map
  before using it, author takes no responsibility. Using the geomag
  library for WMM2010 magnetic declination model.

This would be perfect for what I need. Unfortunately though it only works with QGIS up to 2.18.

Is there anything similar for QGIS 3.X?


texts: Merging DEM file generted many negative value in QGIS

texts:
With four DEM files downloaded by Python Elevation package, I merged them in QGIS 2.8.3.

I use raster-> Miscellanaous -> Input direction, no more attribute was added.

The four dem data has the elevation range roughly from (-100,2000) meters.

But in the merged file, the data range was from -30000 to 2000 meters.

The data distribution are shown here:


texts: Bug in F******** spatial relation? (I'm not cussing)

texts:
I would expect ISpatialFilter.SpatialRelDescription to behave consistently with IRelationalOperator.Relation.  However, when I run the code below, I get:

    Total features : 44368
    Count1: 9 T********
    Count2: 9 T********
    Count1: 21 F********
    Count2: 44359 F********
    Count1: 1 ******F**
    Count2: 1 ******F**

Is this a bug in the F******** relation ?

private void TestSpatialFilter()
{
    IMxDocument mxdoc = ((IApplication)m_application).Document as IMxDocument;
    IFeatureLayer fLayer = mxdoc.FocusMap.get_Layer(0) as IFeatureLayer;
    IEnvelope env = ((IActiveView)mxdoc.FocusMap).Extent;
    env.Expand(0.5, 0.5, true);
    env.Project(((IGeoDataset)fLayer.FeatureClass).SpatialReference);
    Debug.Print("Total features : {0}", fLayer.FeatureClass.FeatureCount(null));
    IFeatureClass fc = fLayer.FeatureClass;
    string[] rels = { "T********", "F********", "******F**" };
    foreach (string rel in rels)
    {
        Debug.Print("Count1: {0} {1}", Count1(fc, env, rel), rel);
        Debug.Print("Count2: {0} {1}", Count2(fc, env, rel), rel);
    }
}

public static int Count1(IFeatureClass fc, IGeometry geom,  string relDescription)
{
    ISpatialFilter sf = new SpatialFilterClass();
    sf.SpatialRel = esriSpatialRelEnum.esriSpatialRelRelation;
    sf.SpatialRelDescription = relDescription;
    sf.Geometry = geom;
    int i = fc.FeatureCount(sf);
    Marshal.FinalReleaseComObject(sf);
    return i;
}
public static int Count2(IFeatureClass fc, IGeometry geom, string relDescription)
{
    int i = 0;
    string relation = String.Format("RELATE(G1,G2,'{0}')", relDescription);
    IRelationalOperator relOp = geom as IRelationalOperator;

    IFeatureCursor fCur = fc.Search(null, false);
    IFeature feat;
    while ((feat = fCur.NextFeature()) != null)
    {
        if(relOp.Relation(feat.Shape,relation))
            i++;
    }
    Marshal.FinalReleaseComObject(fCur);
    return i;
}


texts:
I am new to ArcGIS and encountering the above error code when I'm attempting to extract data from a series of rasters.  Here's some more information:

The process: I overlaid a shapefile containing a number of coordinate points on top of a series of rasters that are tif files.  I didn't modify them because the coordinate system is the same (NAD83).  I did not apply a projection to the shapefile, but the rasters automatically have a projection applied (USA Contiguous Albers Equal Area Conic USGS Version). I'm not quite sure if this is being projected on the fly, because in ArcCatalog, it also has the projection applied.  I then attempted to extract the data from these points using the "Multi values to point" tool.

The issue:  Every time I try to extract the data, I get an error that simply says "Extract multi values to point [000000-12301899]".  I attempted to change the rasters from tif to grid files, and that didn't help either (it executed, but the resulting attribute table had no values).  What am I doing wrong?