texts: ArcGIS Desktop registry value under software key?

texts:
Possible Duplicate:
How can I programmatically get the path of “Python.exe” used by ArcMap

I am trying to get some information from registry about ArcGIS. I am using ArcGIS Desktop 10.0 and For this the registry looks like:

My questions are:

What does it mean when the "Python10.0" subkey has a "True" value?
When "Python10.0" subkey has a "False" value?
Does the registry key values same for every version of ArcGIS?

[Close]
I have edited How can I programmatically get the path of “Python.exe” used by ArcMap question. Thanks @blah238 and @whuber♦ to get me the wrong. Now, this is a close question. Please help me for previous one.

Closely related
At How can I programmatically get the path of "Python.exe" used by ArcMap I asked how to obtain the path to Python.exe programmatically within ArcMap. Of the several solutions offered, one of them suggested using the registry, prompting the present question.


texts: I want to clip out 500 - 1000 cut-outs of 500px*500px of an image (a city area) (.tiff) and a shapefile layer above it (.json or .csv). These files are big for a manual process as I have to select the area manually and clip it. Is there any way this can be automated to cut the image into chunks of 500x500 pixels?


texts:
I've just installed qgis via OSGoe4W custom install on a brand new windows 10 machine.
in the desktop section of the install dialog I've selected grass 7.8.3-1
in libs section QGIS-GRASS-plugin7 for QGIS 3.14.16-3 as well as for QGIS ltr version 3.10.10-3
the installation ran without any troubles.
trying to run any GRASS function via processing toolbox returns the following error:
"the specified GRASS 7 folder  "C:\OSGEO4~1\bin\bin" does not a valid set of GRASS 7 modules. Please, go to the Processing settings dialog, and check that the GRASS 7 folder is correctly configured"
searching StackExchange returned some promising results - like GRASS not working with QGIS 3.7 (and it's predecessor-article for QGIS 3.6 as well) where one user suggested tweaking QGIS settings via "Settings -> Options -> Advanced -> Configuration" and set the current Grass7 path explicitly.
but my problem is, that neither QGIS ltr/3.10 nor QGIS 3.14. do show a "GRASS Folder" entry in the Advanced Settings Editor - see the following screenshot:

so I do not have any option to identify the GRASS path in the same way as described in the before-mentioned posting.
the remarkable thing is that version 3.14. on mac osx looks different:

and ist works perfectly ...


texts: Error installing Photo2Shape QGIS 2.6 for Mac

texts:
I am trying to install the Photo2Shape plugin for QGIS on a Mac (10.9 Mavericks) but I am getting the following error:

The plugin is broken. Python said:
No module named exifread

I have installed exifread using pip install exifread from inside Terminal, but the plugin is still having issues finding the module. Is there a way to link modules inside of QGIS? Or is there other options?


texts: Read a Raster Band Value inside of a VectorLayer in QGIS

texts: How to read pixel values in RasterLayer that intersects line features of VectorLayer, applying Python in QGIS?


texts: Finding intersection length between two polygons using PyQGIS?

texts:
I have two polygons and i want to find the intersection line length between these polygons.

for geom1,geom2 in itertools.permutations(selection,r=2):
    intersection =geom1.geometry().intersection(geom2.geometry())

it gives the boolean values as true/false. But i want to find the intersection length.How can i do that using python?


texts: How to create something similar to GDB in QGIS?

texts:
I used to work with ArcGIS. There I had geodatabases in which I could store several shapefiles with defined layer styles.

Is there a way to do something similar with QGIS?

The background is: we are a team of 5-10 people, and would like to save all of our spatial data in one database, so that we can use them in different QGIS projects on several workstations, and when someone changes something it will be automatically updated in all projects. Before I started here, we had many different "local" shapes, and sometimes changes would effect other peoples' work without them getting to know the changes.

Also it was quite comfortable that in ArcGIS you could just add a geodatabase to a project and then have dozens of shapes inserted at once, with defined layer styles.
At the moment we are saving our shapefiles on a hard disk in our network, and have to add the shapes manually to our projects, but then we have to sort the shapes in different folders and adjust the layer styles manually every time..

Thanks in advance.

Sorry for grammar mistakes, not a native speaker here ;)


texts:
I'd like to edit the nodata values in several raster files, but big number of these file doesn't have any nodata value defined, there is nothing to edit.
I'v tried the gdal_translate -a_nodata command but it'as not working for me.
As I understand this command is only suitable for changing an already defined nodata value.

Any ideas?


texts: Exceptions for invalid files/ filenames that already exist in database

texts:
The script below is functional, but it runs into issues if it comes across duplicate files (even when overwriteOutput=True) and corrupt shapefiles when running feature to feature.

Are there exceptions I can add that will either skip these features if they already exist in the database, or continue running the script if it runs into an error adding a shapefile to a geodatabasing?

This is the script that I am trying to add the error handling to. The script performs the following actions:

Walk through a directory and compile a list of shapefiles
Add the shapefiles to a geodatabase

import arcpy, os, fnmatch
from arcpy import env

env.overwriteOutput=True

rootDir = "whatever directory"
search = "whatever criteria"
outGDB = "whatever GDB"

shpLst = []
    for root, dirs., files in os.walk(rootDir, topdown=True):
        for filename in fnmatch.filter(files, search):
            shpLst.append(os.path.join(root, filename))

for fc in shpLst:
    date = os.path.basename(fc).split("_")[0]
    outName = os.path.basename(fc)[9:].replace(".shp", "_") + date
    arcpy.FeatureClassToFeatureClass_management(fc, outGDB, outName, "", "", "")
    print outName + " sent to " + outGDB


texts: ApplySymbologyFromLayer does not work in my python code

texts:
I have  a layer file. I make a button in ArcGIS 10 with .NET. I want to change layer symbology with the layer file. For this purpose, I write python code and execute it from C#. If the following code is written in python window of ArcGIS, then works fantastic. But from my code, it does not work. Can anyone help me?

try
           {
               const string prg = @"import sys
import arcpy
from arcpy import env
mxd = arcpy.mapping.MapDocument('CURRENT')
layer = arcpy.mapping.ListLayers(mxd)[0]
    arcpy.ApplySymbologyFromLayer_management(layer,'C:\Users\DELL\AppData\Local\Temp\Junctions_Temp.lyr')
"
               ;
           var sw = new StreamWriter("D:\\test2.py");
           sw.Write(prg); // write this program to a file
           sw.Close();
           var p = new Process
                       {
                           StartInfo =
                               {
                                   FileName = @"C:\Python26\ArcGIS10.0\python.exe",
                                   RedirectStandardOutput = true,
                                   UseShellExecute = false,
                                   Arguments = "D:\\test2.py "
                               }
                       }; // create process (i.e., the python program
           p.Start();


texts: Can I use Mapbox to create a web form that finds parks near a location?

texts: I have run a network analysis in ArcGIS and have all of the trimmed polygons along roads within a half a mile of a park.  I need to find a way to create a web form that takes in an address and spits out parks within a half mile of the specified location.  Can I do this using Mapbox?


texts: Making lines from points like Points2One for QGIS 3

texts:
Previously I have been using Points2One to make lines from points in QGIS 2.18 but the plugin is not available in QGIS 3.0.

I am wondering if there is an alternative?


texts: WeedMap pro users - Can weed points and polygons be extracted for use in other GIS clients?

texts:
Are there any users out there of the WeedMap Pro product? I know a lot of local councils in Australia use this for weed management.

I am involved with a regional weeds committee involving a number of councils. Each council is putting forward their weeds infestation layers to produce a region wide map, however the councils are using WeedMap pro and apparently they cannot extract their data for use outside of WeedMap. This doesn't sound right? You would think there has to be some way of extracting the spatial representations of these infestations!

I have been in contact with the WeedMap people (because we are actually looking at the product ourselves - its quite cool, but this kind of issue annoys me!) and they said the product uses a SQL server 2005 and stores the geometry as text, not a db specific spatial data type. Also that they don't allow direct access to the database due to security.

Is there anyone out there using WeedMap pro that knows a way around this? Even if all I end up with is a list of coordinates for points, and polygon centroids that would be acceptable.

Thanks

Ando


texts: QGIS flowmapper migration arrows from origin X to destination X - possible?

texts:
I have a question about QGIS with the Flowmapper plugin. I am trying to visualize some real estate data. I have data about consumers who are actively searching for housing and what I want to visualize is:

In which neighbourhood does the consumer currently live?
In which neighbourhood would the consumer like to find a home?

I have used QGIS to find neighbourhood centroids. After this I created an 80x80 matrix and plotted all the ''migration'' data using the Flowmapper plugin and everything worked out fine.

There is just one thing that I cannot visualize and that is: Consumers who would like to find a home in the same neighbourhood as they currently live in.

The Flowmapper plugin only seems to plot the migration arrows from neighbourhood X to all the other neighbourhoods but not from neighbourhood X to neighbourhood X.

Does anyone know if there is a way to visualize this? (preferably with Flowmapper)

Perhaps with a circular arrow that points from neighbourhood X back to neighbourhood X.

For a complete analysis I would very much like to know how consumers aim to move from one to another neighbourhood but also how many consumers are contempt with the current neighbourhood they live in and would like to find a new housing within this current neighbourhood.


texts: How to Clip a PostGIS raster, and get a GeoTIFF out?

texts:
My end goal is to load the DEM of the entire country in the PostGIS Database, and create a service, which will take in a GeoJSON polygon, and give you a tiff file of that area.

For this, my Research is currently stuck in clipping the PostGIS raster table.

This is what I have done so far:

I loaded a small DEM tiff into the PostGIS database, using raster2pgsql. This has created a Table with 16 rows, and it can be successful seen as a Raster in QGIS.
I checked if the Given polygon intersects with the Table, using the following command:

SELECT ST_Intersects(
    ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57
     ] ] ] }'), 4326),
rast) from carto_dm1;

This showed me that 2 rows of the table intersect with the given geometry, & 14 rows do not.

If I use ST_Clip, I get multiple records back. I'm using the following command:

Select ST_Clip(rast,
    ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57] ] ] }'), 4326)
    ) from carto_dm1
    where
ST_Intersects(
    ST_SETSRID(ST_GeomFromGeoJSON('{ "type": "Polygon", "coordinates": [ [ [ 73.53, 18.57], [ 73.50, 18.45], [ 73.76, 18.48], [ 73.72, 18.58], [ 73.53, 18.57] ] ] }'), 4326),
    rast);

Now how do I get one clipped Raster out, which can be passed to  ST_AsTIFF, so that I can write a GeoTiff from it?


texts: Adding Model Only Tool (Select Data) to model in ArcGIS Desktop?

texts:
I know the Select Data tool exists because of its documentation here:

http://pro.arcgis.com/en/pro-app/tool-reference/modelbuilder-toolbox/select-data.htm

However, when I search for the tool in ArcGIS Desktop 10.5 it does not appear. In fact, I can't find where any Model Only Tools are located.

Is there some trick to finding these like turning them on or downloading them separately?


texts: How to run (SQL QUERY) in a different schema in QGIS DB Manager?

texts:
My last question was How to QUERY (SQL) in QGIS DB Manager?
The result was this:

I was successful...

But when I try to use the same layer in another SCHEMA (Not in Public) it gives me this error:

Could somebody tell me: what i´m doing wrong here?


texts:
Is there a way to extract polygon lines from a map loaded as a .tiff file on QGIS? The original map (see below) of Venice is a scan from a book and has all the lines I need (they mark historical parishes in the city), but it would take to long to manually trace them. I need each parish to be its own entity/polygon.

This image is not georeferenced, and I don't need it to be right now actually. I just find it useful to work on QGIS when dealing with maps.


texts: Using specific ODBC connection when writing to SQL?

texts:
I am currently writing an output file that is generated using some Mapbasic Code into a MS SQL Database:

Dim SQL_output_tablename as String

SQL_output_tablename = """dbo"""&"."&"""ClientOutput_"
   + CurrentDay + "_" + CurrentMonth + "_" + CurrentYear
   + "_" + CurrentHour + "h" + CurrentMinute + "m"""

Commit Table Client As "C:\Program Files\MapInfo\Mapinfo_SQL\Client.TAB"
   **Type ODBC Connection 1**
   Table SQL_output_tablename Type SQLServerSpatial
   Geometry ConvertDateTime ON Interactive

The issue is that I am unsure of how to reference the ODBC connection using the name instead of the assigned Index. Mapinfo gives the connections an Index, 1 being assigned to the first ODBC connection opened etc.

Is it possible to write the code differently and specify the ODBC connection (The name and not the index) to be used when writing to a DB?


texts: In arcmap how to modify time value in editor tracking fields

texts:
In my map, I have several layers enable editor tracking, but the time zone is not the same as our location time zone (GMT +8). When I try to use time slider to modify the time, I find the time slider is disabled.

I want to check how to modify the time value in the editor tracking fields. If use time slide, how can I enable time slider?

I am using arcmap 10.2.2.


texts: Building multiple CON statements with OR/AND in ArcGIS Raster Calculator?

texts:
I have two raster "T" and "N" which have values ranging from 0 to 1000. This is the logic that I want to query and create the new raster:

This is how I tried to query in a syntax-correct way in the Raster Calculator in ArcGIS, but the resulting raster only shows the value "1", but no error.

Con(("T.tif" < 50) & ("N.tif" < 40), 100, 0) |
Con(("T.tif" < 50) | ("N.tif" < 40), 80, 0) |
Con(("T.tif" < 55) & ("N.tif" < 45), 60, 0) |
Con(("T.tif" < 55) | ("N.tif" < 45), 40, 0) |
Con(("T.tif" < 60) & ("N.tif" < 50), 20, 0)

What am I doing wrong?


texts: Unable to send successful WFS request to Geoserver to remote host

texts:
I use Python Requests in order to edit/add spatial data stored in a postgis and served through geoserver.

In the geoserver I have two different stores. One which connects to the server (host is localhost; the server the app runs) and another one with another IP address.

Also the schema is different for each store (one is public).

Everything works fine when I try to send a POST request and edit the location of some point which is stored in the localhost database. But when I try to do the same for points stored in the remote database although I get a success response (200), the edits/additions don't take place.

I was wondering if this is related with the different schema name. Or because the store has a different host. Although that should be handled by geoserver.

Any ideas what can it be the issue?

EDIT

This is part of the logfile from geoserver:

   Request: getServiceInfo
   2016-09-26 17:01:52,588 ERROR [geoserver.ows] -
   org.geoserver.wfs.WFSException: java.lang.StringIndexOutOfBoundsException: String index out of range: -1
   at org.geoserver.wfs.Transaction.transaction(Transaction.java:111)
   at     org.geoserver.wfs.DefaultWebFeatureService.transaction(DefaultWebFeatureService.java:170) ...

And this is the XML request I sent:

<wfs:Transaction service="WFS" version="1.1.0"
 xmlns:ogc="http://www.opengis.net/ogc"
 xmlns:wfs="http://www.opengis.net/wfs"
 xmlns:gml="http://www.opengis.net/gml">
<wfs:Update typeName="geonode:wld_poi_warehouses_wfp_1">
  <wfs:Property>
   <wfs:Name>the_geom</wfs:Name>
    <wfs:Value>
     <gml:Point>
       <gml:pos>48.4583518828 -2.109375</gml:pos>
     </gml:Point>
   </wfs:Value>
 </wfs:Property>
 <ogc:Filter>
   <ogc:FeatureId fid="wld_poi_warehouses_wfp_1.2"/>
  </ogc:Filter>
 </wfs:Update>
</wfs:Transaction>


texts: Categorical raster classification with python in qgis

texts:
I have 5 categorical rasters and an old script to classify them. The script (*.sml) was created for TNTmips which is no longer an option. I know how to do the operation with raster calculator in QGIS but it is a tedious job. The script looks something like this:

    GetlnputRaster(A)
    GetlnputRaster(B)
    GetlnputRaster(C)
    GetlnputRaster(D)
    GetlnputRaster(E)

    GetOutputRaster(F,Numlins(A),NumCols(A),RastType(A))
    SetScale(F,LinScale(A),ColScale(A))
    CopySubobjects(A,F)

    for each A begin
    F=((A=1   and   B<=40 and C<=60)*6\
        +(A=1   and   B<=40  and C>60)*2\
        +(A=1   and  B>40  and  B<=60 and C<=60)*11\
        +(A=1   and  B>40  and  B<=60 and C>60)*21\
        +(A= 1 and  B>60)*1\
    +(A=2  and  C<=60 and B<=40)*61\
    +(A=2  and  C<=60 and B>40)*11\
    +(A=2  and  C>60 and B>40)*11\
    +(A=2  and  C>60  and B<=40)*21\
        +(A=41   and  D<=40 and  B<=50 and  C<=60)*61\
        +(A=41   and  D<=40 and  B<=50 and C>60)*21\
        +(A=41   and  D<=40 and  B>50 and  B<=60 and C>60)*21\
        +(A=41   and  D<=40 and  B>50  and  B<=60 and C<=60)*11\
        +(A=41   and  D<=40 and B>60)*11\
        +(A=41  and  D>40  and  D<=60 and B<=60)*411\
        +(A=41  and  D>40  and  D<=60 and B>60)*11\
        +(A=41   and D>60)*411\
    ..............
    end

but with a total of 54 lines.
Is there a simpler way of doing this in python or to create a text file with classification rules and run it in python or more appropriate, how to translate this script in python?
I ask this because I am not so good in python programming and I want to avoid a very, very long formula in raster calculator.

I have found a possible solution which raise another question.
I have used the Combine tool from ArcGIS for three rasters and then I have added a new field in which I have used a python expression (in Field Calculator) to recode the cells:

def Classify(A, B, C):
  if((A == 100 and B == 1 and C == 1) or
     (A == 100 and B == 2 and C == 1) or
     (A == 100 and B == 5) or
     (A == 200 and B == 1 and C == 1) or
     (A == 200 and B == 2 and C == 1) or
     (A == 400 and B == 5 and C == 1) or
     (A == 500 and B == 5 and C == 1)):
      return 111
  elif ((A == 100 and B == 1 and C == 2) or
        (A == 100 and B == 2 and C == 2) or
        (A == 100 and B == 3 and C == 1) or
        (A == 200 and B == 1 and C == 2) or
        (A == 200 and B == 2 and C == 2) or
        (A == 400 and B == 5 and C > 1) or
        (A == 500 and B == 5 and C > 1)):
        return 112
  elif ((A == 100 and B == 1 and C == 3) or
        (A == 100 and B == 2 and C == 3) or
        (A == 100 and B == 3 and C == 2) or
        (A == 100 and B == 4 and C == 1) or
        (A == 200 and B == 1 and C == 3) or
        (A == 200 and B == 2 and C == 3)):
        return 113
............................................

I have used a slightly different code system but the principle is the same.
Now the question:
There is any possibility to create a table with all the criteria and use this table to recode a raster? I ask this question because I intend in the future to use more than three categorical rasters.


texts: QGIS Vector-/ Raster Import of multiple non-spatially registered datasets

texts:
When you import multiple spatially unregistered datasets and you want to set the same coordinate system to all data you have to set it for every dataset separately in the import dialogue.

Is there a way to remember your decision so that you have to specify the coordinate system only once?


texts: Ignore features outside the area of interest in QGIS for an ArcGIS Feature Service layer

texts:
I have layers in ArcGIS Feature Service format in my map that have more than 100,000 features. I am not editing them but just using them as a reference. I work only in a particular small area at one time. I can't open the attribute table for these files, because they have thousands of features.
Is there a way to set a filter in a QGIS project to ignore everything that is outside a particular area, so that the map is faster, and I can open the attribute table. If the filter is not possible, then only open the attribute table for selected features.


texts: What is the OGC way of representing Related features?

texts:
I have a point shapefile of Persons (The point is located at the registered address of the person) and a polygon shapefile of parcels.

There is an owner_id attribute in the parcels shapefile, which links the owner in the persons shapefile to the the parcel ( For the sake of argument lets assume that one parcel can have only one owner, while one person can own multiple parcels)

If I was working with ArcGIS, I would have created a relationship class to link these two featureclasses in the Geo-database, and then Published the data using ArcGIS Server. I could then use Query Related Records (Operation) on the Rest End point in my Web App.

But I'm working with the Open source stack (Qgis+PostGIS+Geoserver+OpenLayers). What is the way of modelling these relationships in this stack?

Right now I have the primary-key foreign-key relationship set up in the PostGIS Database, and have published these two layers as two different services. Whenever I need to find related records, I have to manually run WFS requests on the other layer.

Is there a better, more standard compliant or best-practice way of doing this in the Open Source Stack?


texts: How to correctly merge layers with overlapping polygons in QGIS/GRASS?

texts:
I have 4 vector shapefile that I would like to join/merge, without worrying about attributes.
Here are the shapefiles.
I would like to do an union and a spatial dissolve, but everytime I do a dissolve after "PUP" and "PRG" layers union, I got an error about wrong geometries and the output layer is empty.
I validated layers geometry with no errors (at least is what QGIS tells me).

Am I missing out something?


texts: ERROR 000157: Input and target dataset should have the same number of bands. Failed to execute (MosaicToNewRaster)

texts:
This is just a continuation of my former post:Batch Mosaic Rasters: arcpy.MosaicToNewRaster_management (ERROR 000628: Cannot set input into parameter coordinate_system_for_the_raster) --where I batch process mosaic rasters using arcpy.MosaicToNewRaster_management. The first file from the "list_all_rasters" has been successfully processed but the next file was not. Pretty weird, since all of my input rasters have the same number of bands and other properties.

This is my code:

   #This is for mosaicing rasters
   import arcpy
   from arcpy import env
   env.overwriteOutput = True

   #Set the current workspace and other env 't variables
   env.workspace = r"C:\thesis\for_sampling\sampling_outputs\new_loss"
   env.nodata = "MINIMUM"
   env.compression = "LZ77"

   list_all_rasters = ["newlossGreater_Luzonforests_onlyloss_20N_120E.tif;newlossGreater_Luzonforests_onlyloss_30N_120E.tif;\
newlossGreater_Luzonforests_onlyloss_20N_110E.tif","newlossGreater_Palawanforests_onlyloss_20N_120E.tif;\
newlossGreater_Palawanforests_onlyloss_10N_110E.tif;\newlossGreater_Palawanforests_onlyloss_20N_110E.tif",\
"newlossGreater_Negros_Panayforests_onlyloss_10N_120E.tif;newlossGreater_Negros_Panayforests_onlyloss_10N_110E.tif",\
"newlossGreater_Mindanaoforests_onlyloss_10N_120E.tif;newlossGreater_Mindanaoforests_onlyloss_10N_110E.tif;\
newlossGreater_Mindanaoforests_onlyloss_20N_110E.tif"]

   output_list = ["Luzon_loss_try.tif","Palawan_loss_try.tif","Negros_Panay_loss_try.tif","Mindanao_loss_try.tif"]

   for raster in range(0,3):
       env.workspace = r"C:\thesis\for_sampling\sampling_outputs\new_loss"
       print raster # checking the list
       arcpy.MosaicToNewRaster_management(list_all_rasters[raster], r"C:\thesis\for_sampling\sampling_outputs", \
                                   output_list[raster], "GCS_WGS_1984.prj", "1_BIT", "0.00025", 1)

   print "Finish all of them!"

And this is the error message:

   ExecuteError: Failed to execute. Parameters are not valid.
   ERROR 000157: Input and target dataset should have the same number of bands
   Failed to execute (MosaicToNewRaster).

I'm pretty sure that that my input files are set to have the same number bands with my output files ( which is 1, in this case). I've double checked my files, and all of them are not corrupted.


texts: Auto-fill attribution of similar feature classes via connectivity rules (Geometric Network)

texts:
I am working with water utility data.  I am trying to figure a way to make feature types connected to the Mains (primary distribution water pipes) automatically fill certain attributes when they are created that already been filled in the Main record.

For example I want to have new service connections (pipes to houses & businesses) that are connected to the Geometric Network (GN) via a Main, to grab certian attribution from the main that the two features have in common.    For instance service area.  And, I want to have this happen when the editor creates the new and connects to the GN.

I am thinking that there should be a way to use the connectivity rules to do this.

Perhaps some sort of python modification to the Transfer Attributes Tool where the search distance is replaced with connectivity.


texts: Printing Soviet-era Russian military maps at original scale

texts:
I have some Soviet era maps of my country that are very detailed. I have access to a plotter and would like to make a hard copy for my hiking trips but have issues getting the resolution/map scale to fit just right.

What app can I use to get the JPEG image as close to its scale as possible?


texts: Conversion of point grid to a raster with open source tools (Gdal, QGIS ...)

texts:
I'm looking to do something which seem simple to me, but I haven't found a solution since 2 days I looked for.

I've got a Postgis DB with lot of points (WGS84) spaced of 100 meters which cover a territory. For each of these points, I've got a value associated. I want to convert all theses point to a raster Geotiff.

The ideal process will make a Geotiff file where each pixels corresponds to the value of the corresponding points of the db.

I tried the heatmap plug-in of QGIS, but It doesn't do what I want, it interpolate lot of new points I don't need, I want in the tiff file only values I've got in my database.

Any suggestion about how to easily generate this file ?


texts: QGIS TIN Interpolation won't run (reports "Wrong or missing parameter value: Input layer(s)"

texts:
I'm trying to run the QGIS Tin Interpolation command, but get the following error message when I try to run the command:

 "Wrong or missing parameter value: Input layer(s)."

The vector input to the command is a shapefile containing less than 2000 points, and a floating point (real) attribute in the shapefile is specified as the interpolation attribute.

Here's a screenshot of the input form:

The error message indicates that the input layer is either wrong or missing:

I've supplied the name of the input vector layer, so the parameter value isn't missing.

Therefore the parameter value must be 'wrong', but it's not clear to me why this would be the case.

Any idea why I'm getting this error message, and what the solution might be?


texts:
Similar question here, I tried the suggested answer but it didn't work for my situation.
updating plugin from qgis 2 to 3. qRegisterResourceData . argument 2 has unexpected type 'str'
I am trying to migrate a plugin to QGIS3 using qgis2to3 and I got stuck on an error;
TypeError: qRegisterResourceData(int, bytes, bytes, bytes): argument 2 has unexpected type 'str'
It seems that when I run the pyrcc5.bat I get a new resources.py but not a new resources.qrc, should have I?


texts: QField sync.FileNotFoundError: [Errno 2]

texts:
I am using QGIS 3.4.5 (Madeira) alongside with QField 1.0.3 (Matterhorn) with the QField Sync 3.1.3.

Everthing worked fine. I made my Project, made it to a QField-Project (with QFiledSync), then to my Tablett, worked in it an then I try to resynchronise it back with my PC. I copied it back from the Tablet to my PC and try to use the Qfield Sync tool to synchronise with my QGIS Project, but every time it gives me the error: FileNotFoundError: [Errno 2] No such file or directory: D:/import Qfield Milan 19/export 8015SW\\data.sqlite.

I never saw that File data.sqlite. What is that for? Why should it be there?  Where is my error?

Even when I pack my Project freshly for QField and directly try to resynchronise the new made QField Project (with no new data) the error occurs.

The full log is:

Ein Fehler trat bei der Ausführung von Python-Code auf:

FileNotFoundError: [Errno 2] No such file or directory: 'D:/import Qfield Milan 19/export 8015SW\\data.sqlite'
Traceback (most recent call last):
  File "C:/Users/Fabia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qfieldsync\dialogs\synchronize_dialog.py", line 65, in start_synchronization
    current_import_file_checksum = import_file_checksum(qfield_folder)
  File "C:/Users/Fabia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins\qfieldsync\utils\file_utils.py", line 89, in import_file_checksum
    with open(path, 'rb') as f:
FileNotFoundError: [Errno 2] No such file or directory: 'D:/import Qfield Milan 19/export 8015SW\\data.sqlite'

Python-Version: 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS-Version: 3.4.5-Madeira Madeira, 89ee6f6e23

Python-Pfad:
C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python
C:/Users/Fabia/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:/Users/Fabia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins
C:/PROGRA~1/QGIS3~1.4/apps/qgis-ltr/./python/plugins
C:\PROGRA~1\QGIS3~1.4\bin\python37.zip
C:\PROGRA~1\QGIS3~1.4\apps\Python37\DLLs
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib
C:\PROGRA~1\QGIS3~1.4\bin
C:\PROGRA~1\QGIS3~1.4\apps\Python37
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\win32\lib
C:\PROGRA~1\QGIS3~1.4\apps\Python37\lib\site-packages\Pythonwin
C:/Users/Fabia/AppData/Roaming/QGIS/QGIS3\profiles\default/python
C:\Users\Fabia\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugins\QuickMultiAttributeEdit3/forms
C:/Users/Fabia/AppData/Local/Temp/tmpwdz6jpz9

I dont know what to do. Do you?


texts: RuntimeError from Python importing ArcPy "The software is not authorized for use"

texts:
In a Windows server there are ArcGIS Pro 2.4 and ArcGIS Server 10.7.1 installed. Both products are licensed.
Moreover, FME server is also installed on that machine.
There are 3 workspaces in the FME server that make calls to Python (ArcGIS server). 2 of them works and the third one is failing.
In order to find the root of the problem, the Python script on the 3rd workspace is run in a command window (for the sake of see it run):

I tried to run the other scripts (pretty similar to this one) and I had the same error.
I open a command window, open Python and run the line:
import arcpy

and got the same result.
After making a quick search, only found this reference to the error:
https://wiki.earthdata.nasa.gov/display/SDT/Run+arcpy+on+ArcGIS+Windows
It suggests that I should use Python from ArcGIS Pro, not from the server. But I must use the one in the server.
The weird, strange, crazy behavior is that 2 of the workspaces work, despite the fact that they also make use of the same Python and both import arcpy too.
Just to have it mind: FME server engines, run under a domain user, and the ArcGIS server runs under other domain user credentials (but I don't think that is related).
Any suggestion?


texts: Turn off the setting of QGIS to produce XML when closing Raster layer

texts:
Every time I open and close a Raster layer in QGIS it creates an XML-file for that raster when closing.
While I do a lot of automatic processing this really bugs me, since I have to add every time when I request a file list to exclude XML-files. I just want QGIS to stop doing this, so do any of you know where to turn this off?


texts: Buffer Until Intersection using ArcGiS ModelBuilder?

texts:
I am trying to use ArcGIS ModelBuilder to create a toolbox that places a buffer of  1 mile around X location, but does not cross major highways. Effectively, I am geo-fencing in real estate comps.

In the model below, for every "active" (property in green) I put a 1 mile buffer to capture all properties "closed" (in blue), then use an intersection to exclude it from crossing the major highways and excluding and closed properties over the boundary.

Unfortunately, this does not work.

Any thoughts on the proper way to do this?

Here is an image of the roads or boundaries that I already have and do not want my buffers to cross.

Here is an image of the current ModelBuilder I have, however it is wrong and does not work.

Here is an image of just a standard 1 mile buffer around X point.

Here is more closely what I am trying to accomplish. Essentially, the logic is buffer up to 1 mile or until you hit a major highway.

Revised ModelBuilder using Smiller's comments

Results of revised ModelBuilder


texts: How to visualize geolocated data on a slippy map on a website?

texts:
I am trying to visualize some data on a website that has been measured while driving around. My source is a CSV file with longitude, latitude and my measured data. I have been able to export it to a GPX file and to display it using OpenLayers using the instructions found on http://wiki.openstreetmap.org/wiki/DE:Openlayers_Track_example.

The problem is that I would like to not only visualize the (points of the) track but also the measured data that has been recorded with them. Some variation in color or line width would be great, but I could not find an example ho to do anything like this in OpenLayers (or another JS library). I guess it would be possible to use rule based styling as explained on http://docs.openlayers.org/library/feature_styling.html#rule-based-styling but I was unable to find out how to retrieve the property from the GPX file.

I found an example that partially matches my needs on geo.dianacht.de/osm/track.html but this uses different data sources for track and data. It would be more useful if I could just use one source. If this should not be possible, any other suggestion is highly appreciated.


texts: Connect to PostGIS using SEXTANTE in QGIS

texts:
I have a localhost PostGIS database that I would like to use with SEXTANTE in QGIS. I can connect to the database just fine using "Add PostGIS layers...", but when I try "PostGIS execute SQL" in SEXTANTE, I´m not sure of what to enter and I can´t find any documentation or examples. All I get is "Wrong database connection name...". Any suggestions/examples on how to connect to the database?

Magnus


texts: Dynamically color labels in color of features (QGIS 3.x)

texts:
How do I dynamically color labels in the color of the labeled feature?
As you can see in the screenshot, the lines are in different colors, but the labels are just in black. As I know, I have to set the label color manually to get my expected result. But is there also a dynamic way to do so? Surely there is, because it‘s QGIS! But I don‘t know how. Any Ideas?


texts: I am trying to install pgRouting 2.0 on postgressql. After writing the data from the share and lib folder into C:\Program Files (x86)\PostgreSQL\9.1... I tried to create pgrouting extension on pgadmin PgSQL console by using the command "CREATE EXTENSION pgrouting". But I got the error that the library librouting.dll doesn't exist. Can anyone help me?


texts: Automaticlally creating polygon shape around multi points of lat longs in shapefile and see it on map?

texts:
I am not a tech guru but we have lat longs of points (of meters). We have created a shapefile, I can see the points on a map.

How do we then automatically draw a shape (polygon) around them all?

Basically creating a map view from a table containing lat / longs.

Assume that I am using a Basic level licence of ArcGIS Desktop.


texts: Is it possible to start the GRASS GUI from R?

texts:
I am using the rgrass7 package in R in order to work with GRASS GIS. After having set the needed enviornment variables, I can now launch a GRASS session from R and and work on my GRASS data using RStudio.

I would like to work with the GRASS GUI and the rgrass7 interface the same way as when I start R (RStudio) from the GRASS terminal (Using RStudio in a GRASS GIS session). The only difference is that I would like to work mainly in RStudio and only open the GRASS GUI only when needed.

Is it possible to launch the GRASS GUI from R within a running GRASS session?

When I try the sollowing command, the GUI does not open and I get a warning message:

> execGRASS("g.gui", ui="wxpython")
Warning message:
running command 'g.gui.exe ui=wxpython' had status 5


texts: Exporting Symbology from ArcGIS for Desktop to Microstation or QGIS?

texts:
I am using ArcGIS 10.1 for Desktop and trying to figure out how to export my symbols of elements in ArcGIS to CAD or QGIS readable format.

The problem is, that I can export symbology to a Style, or export a map to DGN, but there will be an export of just graphics without symbology.

Do you know how to solve this?


texts: Python code to read metadata from raster file

texts: Can you provide me with Python code to read metadata from a raster file?


texts: Definition of custom coordinate system with Equator rotation?

texts:
I need to define cylindrical coordinate system in ArcGIS Desktop with numeric equator rotation. Pictures below show how Equator is rotated.
Coordinates of points 1 and 2 are known.

Is there a way to do it in ArcGIS Desktop?


texts: Are three points really necessary to georeference?

texts:
I find myself in a debate between what I have learned and used, and my employer.

Is it possible to get accurate results georeferencing with only two points instead of three or four?

The answers I have gotten range from:

No.
Yes, if the data is already projected.

I have achieved seamless results georeferencing with two points yet my coworkers think I'm insane or just wrong despite my experience.


texts: Selecting records in table from feature class using ArcGIS ModelBuilder?

texts:
I have a table and a feature class in a file geodatabase. I want to export records from the table that match what is in the feature class using ModelBuilder. For example the table contains 100 records and there are only 30 in the feature class. I want to export only the 30 records that match what is in the table. This is the SQL that I am trying to use:

GIS_ID IN (SELECT GIS_ID FROM FeatureClass)


texts: Possible to use Google Earth to generate "skyline" planes?

texts:
I was really interested to notice that my city appears to have gotten the 3D Buildings treatment on Google Earth. I am trying to figure out how I might these 3D buildings to generate 'skylines' for a given street.

Speaking more abstractly, I think what I want is to generate/export a projection of an arbitrary subset (let's say: 150 > n > 1) of these 3D buildings onto an arbitrary 2D plane of my choosing. More specifically to the goal at hand, one of that 2D plane's dimensions would run normal to the earth's surface, from the surface to some arbitrary height, and the other dimension would run for some arbitrary length of street. (I guess I'll assume for the moment that the street is effectively straight.)

More practically speaking, I envision being able to (in the Google Earth app) place a near clipping plane (i.e. take out all the buildings/terrain/etc between the camera and this plane) and a far clipping plane (i.e. take out all the buildings/terrain/etc behind that plane to infinity), then position the camera and export what the camera sees (again, ideally in a variety of ways, like a standard raster of what I see in Google Earth, or an orthographic wireframe projection, or an orthographic hidden-line projection, or... etc.)

I see that it's possible to download single models, but it seems like it'd be a long haul to get from that to the desired goal, assuming for the moment I even knew what tool I'd use to do the clipping planes/camera/export once I got the models all downloaded. Anyone have any ideas on how I might go about this?


texts: ArcObjects: modifying feature field value

texts:
I'm developing an ArcMap extension and looking for a way to change feature attributes knowing its feature class and OID (and, presumably, everything else that there is in IFeature). All layers use ArcSDE datasets.

What's the right way to do it? Should I change the data in database? It feels like I shouldn't but I can't figure out how to change feature fields values using ArcObjects. The documentation for ArcObjects looks incomplete and awkward to use.

Background info: the extension provides a custom tool that processes selected feature and changes one of its attributes. Now, feature processing involves calling a remote service so it's an asynchronous operation and I should place all my updating logic in the callback method. The tool should work with any of the project layers all of which are ArcSDE datatsets. That's basically the requirements.

If anyone can help me with an advice or a working sample, I would be really grateful. Thanks in advance.


texts: Using specific transformation in R sf::st_transform() to reproduce ESRI’s Project tool

texts:
I am reading in a CSV with latitude/longitude points (WGS 84, EPSG 4326), transforming into the NAD_1983_HARN_StatePlane_Maryland_FIPS_1900_Feet projection (WKID 2893) and saving as a feature class. I wrote an R script to automate work that was being done in ESRI ArcMap.
When performing this operation in ArcMap using the Project tool, the default Geographic Transformation is NAD_1983_HARN_To_WGS_1984_2 (WKID 1900, accuracy 0.10 m).
When I use sf::st_transform(), the results do not match.  R appears to be using the transformation NAD_1983_HARN_To_WGS_1984 (WKID 1580, accuracy 1.0 m). It results in my points being off by ~2.88 feet.
Is there a way to force st_transform() to use the more accurate transformation?
I can't find a way to reference a specific transformation with this package. I even tried building a custom pipeline to pass to st_transform(), but am having trouble understanding the syntax and order of operations.
Here are the State Plane coordinates from my two example data points that show the discrepancy:
R Data:
POINT (1454828 413134.1)
POINT (1454429 485972.8)

ESRI Data:
POINT (1454829 413131.3)
POINT (1454430 485970)

And example code showing how I'm doing the transformation in R
library(sf)
library(arcgisbinding)
arc.check_product()

#create example data and save to CSV
df = data.frame(X_COORD=c(-76.5,-76.5),Y_COORD=c(38.8,39),DESC=c("Pt_A","Pt_B"))
write.csv(df,"df.csv")

#Read CSV as sf object and define the CRS
R_WGS84 <- st_read("df.csv",options=c("X_POSSIBLE_NAMES=X_COORD","Y_POSSIBLE_NAMES=Y_COORD"))
st_crs(R_WGS84) <- 4326

#Transform the data into State Plane projection. This is where things deviate from ArcMap.
R_StatePlane <- st_transform(R_WGS84,2893)

#Write feature class to file GDB
arc.write("PathToExampleGDB.gdb/R_StatePlane"), R_StatePlane,shape_info=list(type='Point',WKID=2893),overwrite = TRUE)

And an example of the "Project" Tool in ArcMap


texts: openlayers + mapserver: controlling tiles

texts:
I have a mapserver serving the following mapfile:

MAP
    NAME "GLOBCOVER MAP"
    # Map image size
    SIZE 1000 1000
    UNITS meters

    EXTENT -30.001389 -30.106069 80.998611 80.108847

    PROJECTION
        'proj=longlat'
        'ellps=WGS84'
        'towgs84=0,0,0,0,0,0,0'
        'no_defs'
    END

    # Background color for the map canvas -- change as desired
    IMAGECOLOR 255 255 255
    IMAGEQUALITY 95
    IMAGETYPE agg

    OUTPUTFORMAT
        NAME agg
        DRIVER AGG/PNG
        IMAGEMODE RGB
    END

    # Legend
    LEGEND
        IMAGECOLOR 255 255 255
        STATUS ON
        KEYSIZE 18 12
        LABEL
        TYPE BITMAP
        SIZE MEDIUM
        COLOR 0 0 89
    END
END

# Web interface definition. Only the template parameter
# is required to display a map. See MapServer documentation
WEB
    # Set IMAGEPATH to the path where MapServer should
    # write its output.
    IMAGEPATH '/tmp/'

    # Set IMAGEURL to the url that points to IMAGEPATH
    # as defined in your web server configuration
    IMAGEURL '/tmp/'

    # WMS server settings
    METADATA
        'ows_title'           'QGIS-MAP'
      'ows_onlineresource'  'http://localhost/cgi-bin/mapserv?map=/opt/local/apache2/htdocs/eth_luc.map'
  'ows_srs'             'EPSG:4326'
END

#Scale range at which web interface will operate
# Template and header/footer settings
# Only the template parameter is required to display a map. See MapServer documentation
TEMPLATE 'fooOnlyForWMSGetFeatureInfo'
END

LAYER
  NAME 'GLOBCOVER_L4_200901_200912_V2.3'
  TYPE RASTER
  DUMP true
  TEMPLATE fooOnlyForWMSGetFeatureInfo
EXTENT -30.001389 -30.106069 80.998611 80.108847
  DATA '/Users/calvin/work/luc/proj_public/media/GLOBCOVER_L4_200901_200912_V2.3.tif'
  METADATA
      'ows_title' 'GLOBCOVER_L4_200901_200912_V2.3'
  END
  STATUS OFF
  TRANSPARENCY 100
  PROJECTION
      'proj=longlat'
      'ellps=WGS84'
      'towgs84=0,0,0,0,0,0,0'
    'no_defs'
   END
END

END

A png is dynamically generated by mapserver and shown on the web browser as expected from what is specified in my mapfile's EXTENT.  The url can be seen in the screenshot here:-

However, when I attempt to make use of openlayers to do the same using the same mapfile, by means of this openlayers code snippet:

<script type="text/javascript">
    $(document).ready(function() {
        var lon = 5;
        var lat = 40;
        var zoom = 5;
        var map, layer;
        var defaultProjection = new OpenLayers.Projection("EPSG:4326")

        map = new OpenLayers.Map({
            div: 'map',
            displayProjection: defaultProjection
        });
        layer = new OpenLayers.Layer.MapServer(
            "GLOBCOVER_L4_200901_200912_V2.3",
            "http://localhost/cgi-bin/mapserv",
            {
                map: '/opt/local/apache2/htdocs/eth_luc.map',
                srs: 'EPSG:4326',
                isBaseLayer: true
            });
        console.log(layer);
        map.addLayer(layer);
        map.setCenter(new OpenLayers.LonLat(lon, lat), zoom);
        map.addControl( new OpenLayers.Control.LayerSwitcher() );
        console.log(map)
    });
</script>

all I am getting are 36 instances of white PNGs tiled and displayed on the browser, as seen here:-

What am I missing in order to serve my dynamically generated mapserver tiles (PNGs) on openlayers? (with the same visual result as what was done by the pure mapserver solution)


texts: Making rivers with higher Strahler/stream order wider in QGIS?

texts:
I have a river network and the Strahler order of all rivers in it. Now I would like to create a map similar to this one http://www.hydrosheds.org/, i.e. with rivers of a higher Strahler order being wider than those with a lower Strahler order.

How can I perform that kind of analyses?


texts: How to filter fields values collection of selection layer in SQLQueryDialog

texts:
I'm trying to show SQLQueryDialog in ArcObjects for selection layer (layer created from feature layer via layer's context menu Selection -> Create Layer From Selected Features). When query builder is opened the list of fields values consists of all values of the corresponding fields of the layer used for creating selection layer. But in the Select by Attributes dialog in ArcMap values are filtered for selection layers (only values of features of the selection layer are shown).

My question is how to filter fields values in SQLQueryBuilder for selection layer


texts: Assigning random points within polygon using ArcGIS Desktop?

texts: I have a vector layer with multiple polygons. I converted the layer to a raster layer in order to convert the raster cells to points (I needed my points to be in rows). Is there a way to combine the points and polygons, in a way that I could assign all points within the same polygon one value and other points within a different polygon different values for the points. I would like the points within the polygon to have the same values, but different polygons will have different point values.


texts: Finding distance between paired points in separate shapefiles (paired by attribute data)?

texts:
I have a point layer with site_code attribute data. I have another layer of points in different locations but with the same site_code data.

How can I find the linear distance between points that are paired by their site_code attribute data?


texts: Importing large OpenStreetMap database into QGIS?

texts:
I've got a large .osm file I'd like to load into QGIS and stylise with some QML files I downloaded from the web (https://github.com/anitagraser/QGIS-resources/blob/master/qgis2/osm_spatialite).

The file itself is too large to directly import by the add vector layer tool (it's 6gb!), and I find the import from xml/export to spatial lite method quite time consuming, having to check all the tags that I'd like, and doing this process for points, lines and polygons. Perhaps I'm not doing this correctly?

In any case, I have this large .osm file, and I'd like to access the data and stylise based on the QML files I downloaded, what would be a method for doing this?

I'm using a mac.

Edit: it sounds like I need to thin the database out first. I'll play around a little more, and try to do the same on a smaller database and see if I have any luck!


texts:
Let us say that I would like to store latitudes and longitudes like this:

lat
35 degree

long
80 degree

Which system of units does the unit 'degree' belong to? It is not really a SI unit nor a Imperial unit ...


texts:
I have a map service that has more than 60 data layers. For faster performance, I have separated the data with daily updates as a dynamic service. I still have a map service that base data (50+ layers). I have implemented scale rendering to make sure that only relevant content is shown at a given scale. I am still not getting good speed and trying to figure out the best way to make this service even faster.

I am looking for some MAJOR strategies that would allow to make it faster. These services are intended to be used within a Silverlight application.


texts:
I have two raster files:

Raster 1 has 7974 grids (coarser resolution) with Rowid, Value and Count.

Raster 2 is finer resolution and each grid has values either 1 or 0.

What I want to do is to find the number of 1s in each coarser grid (Raster 1). I assumed I could do this using zonal statistics (zonal sum). I used the following settings:

Input raster or feature zone data:
Raster 1

Zone field:
Value

Input value raster:
Raster 2

Statistics type:
Sum

This should just sum up all the 1s falling under the coarser resolution and give me what I need. However, when I run it, it returns me a raster file with only two values: 1 and 0 for each of my coarser grid. What is it that I am doing wrong here?


texts: Making stops not in certain sequence in order to find shortest route with ArcGIS Network Analyst?

texts:
Is it possible to make stops not in a certain sequence in order to find the shortest route with the network analyst?

I now use the small green points as an input stops because I want the lines with a pink dot beside it in the route. Then I reorder them and don't preserve the end and startpoint. I thought this meant that the network analyst doesn't make a real sequence of the input stops and doesn't choose a real start and end point, and just finds the shortest route in between the points.

However if I look at the output of the tool (light green), it seems like the network analyst tool still uses a kind of sequence and begin and endpoint. Because in the picture it looks like the point where the blue arrow points to is a start or end point. The red line is the route I want. So the route makes a round now (at the left where the black arrow is) going past the south segment instead of going to the north  and connect via the small grey link near the blue arrow (like the red route I want). The green output I now receive is not the shortest route.

I thought that one problem could be that the point near the blue arrow is assigned by the tool as start or end point and the other problem that the network analyst takes segment near the black arrow because the point in the south west (black circle) comes later in the sequence than the point in the middle (black circle). So, it gives each point 1,2,3 and so on, instead of tread them like all 0,0,0. So, I thought that this could be solved by having the stops as equal points (all zeros or something) and calculating the shortest route in between those instead of having the stops as a sequence and going from 1 to 2 and 3 and so on. But I don't know if this is right, see also last comment of @ChrisW, Using line feature as stops input in network analyst?.


texts: Removing lines on coast side using ArcGIS Desktop?

texts:
I have a shapefile of European countries (with regional subdivisions) as well as a shapefile of lines.

How do I delete all the lines on the coast-side, using ArcGIS Desktop?


texts:
I am getting sql spatial results from my database query that consists of long/lat polygons.

I am using these results in my c# windows form application.

How can I draw each polygon from their long/lat points in my windows form application?

I want each drawn polygon to be selectable, so that when a user clicks on the shape an event will be raised to perform some task in the code behind.

I used an open source library dotspatial to add a map control on my windows form and adding a world shapefile to show the location of the polygons that will be drawn, but I can't find a way that how can I add a layer to the dotspatial map control to draw polygons.

Is there any other way that I can accomplish this task?

This is a windows desktop standalone application with no internet connectivity.


texts: Merging polygons and having them appear in specific drawing order in ArcMap

texts: This is in ArcMap 10.6.1. I have four polygons that I would like to merge into a single shapefile and have the order of polygons shown on the map from smallest to largest (see first screenshot attached). However, when I try to use the Merge tool, I get a different result (see second screenshot). I'm thinking the merge is done by how it looks at the attribute table of each polygon. Each polygon has the same fields in their attribute tables (see third screenshot for an example). I'm unsure how I can make the Merge tool (or possibly the Append tool) work to order the polygons from smallest to largest like in screenshot 1.


texts:
I try to make a candy cane style for a line geometry with QGIS
I used the hatching and line symbology but the hach are squared.
Have you one better idea to do:

a more smooth and rounded symbology
the hatches limited by the line itself


texts:
I tried to fanout data but it is all in the same FGDB. So, is it possible to do this? I need to write data for every time different amount of different area into different, but schematic same filegeodatabse.

Best regards.


texts: QGIS vector line layer doesn't accept Cyrillic letters

texts:
In QGIS 2.6 vector line layer doesn't accept Cyrillic letters, how can I fix this?

For polygon it it works fine but for line and point it only gives me set of question marks - ?????

It's .shp, I don't know about encoding


texts:
I have a shapefile of road network. I have found out the nodes on the network.

I want to calculate all the possible routes from a start point to the end point. I am using PostgreSQL, PostGIS and pgRouting.


texts:
I have three maps (data frames) in my QGIS composer and would like to link the size, rotation and extents of the maps so that when I 'move item content' or change the rotation in one of the maps it will change in the other maps.

Is there a way to do this?


texts:
I'm trying to do a quite simple operation, an intersection, but when I run it it says : "Input layer A contains invalid geometries (feature 1745). Unable to complete intersection algorithm."

I'm new to GIS so I have little idea what this could mean and where I could search.. How can I determine what is wrong?


texts: What sample size should be used when applying the National Standard Spatial Data Accuracy (NSSDA)?

texts:
I'm less interested in what is recommended by FGDC, and more interested in information like this: SAMPLE SIZE AND CONFIDENCE WHEN APPLYING THE NSSDA, which suggests you need a sample of 100.  Is there any agreement on this in the GIS Community?

I'm not sure if this means anything, but the actual worksheet provided by the FGDC has 40 rows.

Please think twice before just tossing a number into an answer without an explanation. I expect several good answers, but this question isn't intended to be a poll.


texts: How to return FeatureLayer from gp service(using custom function tool) (.NET)?

texts:
My basic problem is that I have a custom function tool HLZCalulation that I created. It takes in a string to a file gdb , opens it, does some calculations, and then returns a featurelayer with the results of said calculations. The tool works as both a simple tool and a model, returning the correct features and not reporting any errors. However, when I publish the model as a service and try to run it from ArcMap, it reports success message from the server but the output feature class is empty.  When I check the server job folder, it has created the shapefile with the features inside it but fails to return the featurelayer without throwing any specific errors (the REST portal returns success message, then a fail message – see error message below).

Packing the parameter is done at the end of execute (once the shapefile is created) rather than validate as with the calculate area sample on the esri site, which could be a problem. The output parameter definition is included below, as well as the lines of code for packing the parameter, and the output from running the rest portal(in browser)

any ideas?

//////////////////////////////////////////////////////////PARAM DEFINITION

                // Output parameter (Derived) and data type is GPFeatureLayerTypeClass
                IGPParameterEdit2 outputParameter = new GPParameterClass();
                outputParameter.DataType = new GPFeatureLayerTypeClass();
                // Value object is GPFeatureLayerClass
                outputParameter.Value = new GPFeatureLayerClass();
                // Set output parameter properties
                outputParameter.Direction = ESRI.ArcGIS.Geoprocessing.esriGPParameterDirection.esriGPParameterDirectionOutput;
                    outputParameter.DisplayName = "HLZ Output FeatureClass";
                    outputParameter.Name = "hlz_featureclass";
                    outputParameter.ParameterType = ESRI.ArcGIS.Geoprocessing.esriGPParameterType.esriGPParameterTypeDerived;
                    parameters.Add(outputParameter);

//////////////////////////PACKAGE PARAM

        IGPValue resultFeatureLayer = m_gpUtilities.CreateFeatureLayer(resultFeaturePath + ".shp");
        IGPParameter parameter = (IGPParameter)userParamValues.get_Element(OUTPUT_PARAM);
        m_gpUtilities.PackGPValue(resultFeatureLayer, parameter);

//////////////Rest Portal output

esriJobMessageTypeInformative: Submitted.

esriJobMessageTypeInformative: Executing...

esriJobMessageTypeInformative: Executing (Model): Model # D:\data\test\TGD_CLIP.gdb

esriJobMessageTypeInformative: Start Time: Mon Oct 11 17:29:18 2010

esriJobMessageTypeInformative: Executing (HLZ Calculation): HLZCalculation D:\data\test\TGD_CLIP.gdb #

esriJobMessageTypeInformative: Start Time: Mon Oct 11 17:29:18 2010

esriJobMessageTypeInformative: Intersect completed, should return

esriJobMessageTypeInformative: Executed (HLZ Calculation) successfully.

esriJobMessageTypeInformative: End Time: Mon Oct 11 17:32:08 2010 (Elapsed Time: 2 minutes 50 seconds)

esriJobMessageTypeInformative: Executed (Model) successfully.

esriJobMessageTypeInformative: End Time: Mon Oct 11 17:32:08 2010 (Elapsed Time: 2 minutes 50 seconds)

esriJobMessageTypeError: Failed.


texts: Using SVG icon as marker for a geoJSON point file in Leaflet?

texts: I have a basemap that displays locations from a geoJSON file. Instead of the default markers I would like to use my own svg files, but can't find out how.


texts:
I have images from UAV which are geo-tagged. I wish to geo-reference images which were shot in sequence, and want to create a mosaic, which further will be a small map.

Is there any possible way in QGIS to do so?


texts: Looking for an alternative to Minimum Bounding Geometry (convex hull) tool for creating polygons from points

texts:
I have lines of points that I'd like to connect to create individual polygons.
The Minimum Bounding Geometry tool is great because it lets me input many points at once and use the "Group" option to create separate polygons. However, the Convex Hull option isn't accurate enough in drawing some of these boundaries. See:

I am looking for an alternative way to turn many points into polygons. Data looks something like this:

Ax1   Ay1   A
Ax2   Ay2   A
...   ...   ...
Bx1   By1   B
...   ...   ...

These points are listed in a particular order, so a Point --> Line --> Feature task would work. I am just looking for an easier way.


texts: Rotating legend item labels in QGIS Layout Manger?

texts: I want to create horizontal legend in QGIS (v 3.8.0). Is it possible to rotate legend item labels in Layout Manager?


texts: How to clean data in attribute table in QGIS 2.18.14?

texts:
I have a 'name' field with composite name; numbers and words
e.g. "1237, Café"

I managed extract the part number with Field Calculator and make a new field (see 1) with this expression

  left(  "Name" , strpos(  "Name"  ,'\s') - 1)

But I can´t find the way to extract the part names - just words!

If I use this

  regexp_replace(  "Name" , '\W', ' ')

Creating a new field result only fields with numbers...(see 2)

EDITED: correction; the results with only numbers; if the type is "whole number (integer)"!!

Updating the 'name' field, result the same fields only without ','... e.g. "1234 Café"


texts: Calculate ellipsoidal area for a projected layer in PyQGIS 3

texts:
I'd like to loop though each feature in a layer and calculate its ellipsoidal area using QgsDistanceArea() in PyQGIS 3.

Everything works fine on unprojected layers, but I get strange results when input layer is in projected CRS.

So far I have written this:

from qgis.core import *

def calculate_area(in_lyr_name, ellipsoid, units):
    input_layer = QgsProject.instance().mapLayersByName(input_lyr_name)[0]

    area = QgsDistanceArea()
    area.setEllipsoid(ellipsoid)

    features = input_layer.getFeatures()
    feature_area = 0
    error_features_counter = 0

    for i, feat in enumerate(features):
        geom = feat.geometry()
        polygon_area = 0
        try:
            if geom.isMultipart():
                polygons = geom.asMultiPolygon()
                for polygon in polygons:
                    polygon_area += area.measurePolygon(polygon[0])
            else:
                polygon = geom.asPolygon()
                polygon_area = area.measurePolygon(polygon[0])

            feature_area += polygon_area

        except AttributeError:  # catches NoneType geometries (broken etc.)
            error_features_counter += 1
            pass

    # calculated area is in sq. metres (see the "else" case)
    if units == 'km²':
        final_area = feature_area / 1e6
    elif units == 'Ha':
        final_area = feature_area / 10000
    else:
        final_area = feature_area
    return final area, error_features_counter

How do I get a constant result that does not depend on the presence of projection?


texts: Export geodatabase relationship to shapefile

texts:
I have a personal geodatabase with a polygon feature class with a one-to-many relationship to a table.

What I want to do is to export to a shapefile with all the infornation associated in every polygon.

I use ArcGis 10.2.1


texts: "OGR: not enough data error" when using st_transform in R

texts:
I'm trying to use st_transform (from the R package sf) on a sf object and am getting the following error when projecting from WGS 1984 to Albers Equal Area:

> sf_obj_prj <- sf_obj %>%
   st_transform(crs =
     "+proj=aea +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs")
OGR: Not enough data
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
GDAL Error 1: Length of input WKB is too small

I have checked that my map has a CRS:

> st_crs(sf_obj)
Coordinate Reference System:
EPSG: 4326
proj4string: "+proj=longlat +datum=WGS84 +no_defs"

I have also checked that all polygons are valid:

> which(!is_true(st_is_valid(sf_obj)))
integer(0)

(is_true is my user-defined function that returns TRUE if true and FALSE if NA/false)

I've tried searching but could not find anything about this error.

EDIT 1:
It seems that this is arising after using rbind to bind two objects.

Here is the st_bbox on the two objects:

> st_bbox(sf1)
      xmin       ymin       xmax       ymax
-73.991244 -11.145219 -66.619911  -7.111824
> st_bbox(sf2)
     xmin      ymin      xmax      ymax
-62.68262  80.03352 -57.99074  82.14633

st_transform works on the objects separately, but not when I use rbind:

> x = st_transform(sf1, crs = prj_alber)
> y = st_transform(sf2, crs = prj_alber)
> z = st_transform(rbind(sf1, sf2), crs = prj_alber)
OGR: Not enough data
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
  GDAL Error 1: Length of input WKB is too small

EDIT 2:
Here is code to reproduce the error (note: requires downloading data from the Brazilian statistics office)

download.file(
  "http://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2000/setor_rural/projecao_geografica/censo_2000/e1000_generico_dxf/uf/ac.zip",
  paste0(directory, "sf1.zip"))
unzip(paste0(directory, "sf1.zip"), exdir = paste0(directory, "sf1"))

download.file(
  "http://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_de_setores_censitarios__divisoes_intramunicipais/censo_2000/setor_urbano/ac/1200203/1200203.zip",
  paste0(directory, "sf2.zip"))
unzip(paste0(directory, "sf2.zip"), exdir = paste0(directory, "sf2"))

sf1 <- st_read(paste0(directory, "sf1/AC/12SE1000G.dxf"))
# (original crs undefined)
st_crs(sf1) <- "+init=epsg:5527"
sf1 <- st_transform(sf1, crs = "+init=epsg:4326")
sf1 <- sf1 %>% dplyr::select(ID = Text)

sf2 <- st_read(paste0(directory, "sf2/1200203.SHP"))
# (original crs was defined for northern hemisphere, this is the epsg code for southern hemisphere)
st_crs(sf2) <- "+init=epsg:32718"
sf2 <- sf2 %>% st_transform(sf2, crs = "+init=epsg:4326")
# (some polygons were invalid so fixing them here)
invalid = which(!is_true(st_is_valid(sf2)))
sf2[invalid, ] = st_make_valid(sf2[invalid,])
sf2 <- sf2 %>% dplyr::select(ID = ID_)

> x = st_transform(sf1, crs = prj_alber)
> y = st_transform(sf2, crs = prj_alber)
> z = st_transform(rbind(sf1, sf2), crs = prj_alber)
OGR: Not enough data
Error in CPL_transform(x, crs$proj4string) : OGR error
In addition: Warning message:
In CPL_transform(x, crs$proj4string) :
  GDAL Error 1: Length of input WKB is too small

But as per Spacedman's solution below, removing the extra dimension fixes the issue:

z <- st_transform(st_zm(rbind(sf1, sf2)), crs = prj_alber)


texts: Force pyramid zoom level on any WMS, WMTS or XYZ?

texts:
Is there a way to lock the view to any tileset level regardless of the scale, e.g remain at zoom level 12 whether the scale is 1:250000 or 1:62500, so that features like text and paths appear at different size rather than adapting to the view ?

I'm on QGIS 3.10


texts: Where can I find data related to hydro power production in india?

texts: This could possibly include meteorological data/precipitation data aggregated on the district level and power line/ electric transmission infrastructure locations


texts:
I built a python addin for Arcmap and I need to print a folder path to a combo box, which should be a destination folder to copy my *.gdb file.

After a click on a "Backup" button this tool should copy my file to a folder on the right. But still I don't know how to figure that out. The only thing I can see is to read the last thing in Arcmap's python window and probably set it as input for a button.

But still have no idea even how to perform it.

EDIT:

I was changing only lines in python script, config.xml was made by default.
Here is a piece of code:

import arcpy
import pythonaddins

class Backup(object):
    """Implementation for BackupServer_addin.Backup (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        if "" in BackupServer_addin.Folder.self.items:
            print self.items

class Folder(object):
    """Implementation for BackupServer_addin.Folder (ComboBox)"""
    def __init__(self):
        self.items = [""]
        self.editable = True
        self.enabled = True
        self.dropdownWidth = 'WWWWWW'
        self.width = 'WWWWWW'
    def onSelChange(self, selection):
        print "New Selection:",selection
    def onEditChange(self, text):
        print "Edit Change:",text
    def onFocus(self, focused):
        pass
    def onEnter(self):
        print "Current Value:",self.value
    def refresh(self):
        pass

This string is just a try to make a reference to the Folder combobox  if "" in BackupServer_addin.Folder.self.items:

Here are the contents of config.xml:

<ESRI.Configuration xmlns="http://schemas.esri.com/Desktop/AddIns" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"><Name>Python Addin</Name><AddInID>{413e6443-7e28-43e0-832a-33730a790ca4}</AddInID><Description>New Addin</Description><Version>0.1</Version><Image /><Author>Untitled</Author><Company>Untitled</Company><Date>02/03/2017</Date><Targets><Target name="Desktop" version="10.2" /></Targets><AddIn language="PYTHON" library="BackupServer_addin.py" namespace="BackupServer_addin"><ArcMap>
    <Commands>
        <Button caption="Backup" category="Python Addin" class="Backup" id="BackupServer_addin.Backup" image="Images\backup.png" message="Backup" tip="Backup"><Help heading="Backup">Backup</Help></Button>
        <ComboBox caption="Folder" category="Python Addin" class="Folder" id="BackupServer_addin.Folder" itemSizeString="WWWWWW" message="Folder" rows="4" sizeString="WWWWWW" tip="Folder"><Help heading="Folder">Folder</Help></ComboBox>
        </Commands>
    <Extensions>
        </Extensions>
    <Toolbars>
        <Toolbar caption="Toolbar" category="Python Addin" id="BackupServer_addin.toolbar" showInitially="true"><Items><Button refID="BackupServer_addin.Backup" /><ComboBox refID="BackupServer_addin.Folder" /></Items></Toolbar>
        </Toolbars>
    <Menus>
        </Menus>
    </ArcMap></AddIn></ESRI.Configuration>


texts: Why my python script of r.watershed has no output？

texts:
I type those code in pyscripter,

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'acc1', drainge = 'dra1', basin = 'bas1', stream = 'str1', half_basin = 'hbas1', visual = 'vis1')
print g.read_command('g.list', _type = 'rast')

PS. I set the output half.basin as half_basin because the following error will occur as I set it as former one,

Message File Name   Line    Position
SyntaxError
    keyword can't be an expression (module2.py, line 27)    C:\Users\Heinz\Documents\pyscript\module2.py    27

and here's the output,

there are no new raster file generated, and I don't know where I make mistake.

Update#1
This is my code:

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'TC'
gsetup.init(gisbase, gisdb, location, mapset)

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'accu1', drainge = 'dran1', basin = 'basi1', stream = 'stre1', half.basin = 'hbasi1', visual = 'visu1')

#check result
print g.read_command('g.list', _type = 'rast')

I use drainage and half.basin as my parameter, and here's the output:

Update#2
I tried the script below,

import os
import sys

#set up GRASS environment variables
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
import grass.script as g
import grass.script.setup as gsetup
gisbase = os.environ['GISBASE']
gisdb = 'C:\Users\Heinz\Documents\grassdata'
location = 'newLocation'
mapset = 'water_test_2'
gsetup.init(gisbase, gisdb, location, mapset)

#import raster data
#g.run_command('r.in.gdal', input = 'hc', out = 'hc')

#run r.watershed module
#set thresh = 5000
g.run_command('r.watershed', elevation = 'hc', threshold = 5000, accumulation = 'accu', drainage = 'dran', basin = 'basi', stream = 'stre', half = 'hbas', visual = 'visu')

#check result
print g.read_command('g.list', _type = 'rast')

And the output looks fine,

*** Remote Interpreter Reinitialized  ***
>>>
----------------------------------------------
raster files available in mapset <water_test_2>:
accu basi dran hbas stre visu

raster files available in mapset <PERMANENT>:
hc

----------------------------------------------

My conclusion is, we can use python to call r.watershed in GRASS 6.4.3 from outside, but one of the input name must be 'half'

Thanks everyone that helps me successfully get the result I want.


texts: OpenLayers 5: requesting WMS layer without CRS like in OpenLayers 2

texts:
Back in the day of OpenLayers 2 displaying a WMS layer involved something like:

var map = new OpenLayers.Map(
    'map',
    {maxExtent: new OpenLayers.Bounds(
         -20037500.000000,-13347182.804674,
          20037500.000000,13347182.804674),
     maxResolution: 66903.171953}
);
var mslayer = new OpenLayers.Layer.MapServer(
    "MapServer Layer",
    "http://localhost/cgi-bin/mapserv?map=/mapfiles/countries.map&",
    {layers: 'Countries'},
    {singleTile: "true", ratio:1}
);
map.addLayer(mslayer);
map.zoomToMaxExtent();

No explicit CRS was declared, and thus the server would respond with the default CRS of the layer.

In OpenLayers 5 a View object is required, in which a projection property is available. In case this property is not set up, all requests made to server explicitly declare the CRS identified with "EPSG:3857" (whether the service supports it or not). To avoid this a CRS must be explicitly declared, registered with proj4js and then passed to the View object:

var Map                 = ol.Map;        //~ import Map from 'ol/Map.js';
var View                = ol.View;       //~ import View from 'ol/View.js';
var TileLayer           = ol.layer.Tile; //~ import TileLayer from 'ol/layer/Tile.js';
var OSM                 = ol.source.OSM; //~ import OSM from 'ol/source/OSM.js';
var TileWMS           = ol.source.TileWMS; //~ import TileWMS from 'ol/source/TileWMS.js';
var proj4             = proj4;             //~ import proj4 from 'proj4';
var register            = ol.proj.proj4.register; //~ import {register} from 'ol/proj/proj4';
      var getProjection       = ol.proj.get;     //~ import {get as getProjection} from 'ol/proj';

proj4.defs('EPSG:152160',
        '+proj=igh +lat_0=0 +lon_0=0 +datum=WGS84 +units=m +no_defs');
register(proj4);
const homolosine = getProjection('EPSG:152160');

var layers = [
    new TileLayer({
      extent: [-1.99496e+07, -8.68326e+06, 1.99675e+07, 8.41393e+06],
      source: new TileWMS({
        url: 'http://localhost/cgi-bin/mapserv?map=/mapfiles/countries.map',
        params: {'LAYERS': 'Countries', 'TILED': true},
        serverType: 'mapserver',
        // Countries have transparency, so do not fade tiles:
        transition: 0
      })
    })
  ];
var map = new Map({
    layers: layers,
    target: 'map',
    view: new View({
      projection: homolosine,
      center: [0, 18000],
      zoom: 4
    })
  });

The problem is that proj4js does not support the same CRSs as PROJ or GDAL. It is easy to arrive at a situation where a layer provided by a server like Mapserver in a certain CRS can not be dealt with in OpenLayers 5.

I would like to understand whether this situation can be circumvented in OpenLayers 5 (the tutorials do not go further than the example above). Hypothetically, this could be achieved in two ways:

a) suppress the explicit CRS declaration and let the server reply with the default CRS of the layer;

b) declare the request CRS using a simple string, avoiding a proj4js registered CRS.


texts: Vector layers seem not to be aligned when performing 'union' in qgis

texts:
I have two files (shapefiles). A file with 50x50 meter cells. and a grid with 100 by 100 meter cells. I want to do 'union' so that I know which 4 cells of the 50x50 meter vector layer make up one cell of the 100 by 100 meter layer.
However, when I do this, there is for some reason a mistake when doing this. I am expecting 4 different ID values for the 50x50 meter file to 1 ID value for the 100x100 meter grid. However, there is many more. This makes me conclude the layers are not aligned. However,

I use the same projection for both vector layers, so I think that should not be the problem.
When I zoom in as far as I can and I change the transparency of the top layer, I still do not see even the slightest mistake in terms of alignment.
I have also tried intersect.


texts: What are the alternatives to ArcMap's Report Wizard?

texts:
I have a polygon layer which has lots of attributes, and I'd like to create a report template which summarises these attributes for individual polygons.

That is, the user will select a polygon from one layer, and see a formatted report showing the relevant attributes for that polygon. There are also features in multiple other featureclasses which are related via an ArcMap relationship, and which should also be summarised on the same report.

It's been a few years since I last used ArcMap's Report Wizard, and it doesn't seem that much has changed in the meanwhile. It still looks quite complicated and cumbersome to set up and format a report, so I'm looking for any alternatives which might exist.

Before I write something from scratch, can anyone recommend any other approaches or code samples which might help? Since the data is stored in a file geodatabase, with the relationships defined in ArcMap, I am assuming I'll need to use arcpy to obtain a reference to the data and traverse the relationships.

StackOverflow has a question on creating an HTML report using Python which will probably be my starting point, unless there's anything else out there.


texts:
I'm trying to merge two vector layer as part of a graphic model but get an error:

  OBJECTID field in layer layer1 has different data type than in other
  layers (Integer64 instead of Real) Execution failed after 0.06 seconds

And both layers have a field called OBJECTID with different data types. I don't care about the attributes I just need the geometries. Isn't there a way of telling merge to ignore the attributes?

I would rather not have to drop/rename all fields prior to merge (I have a lot of layers, not only these two with many fields to drop). Maybe there is a QGIS equivalent of ArcGIS MakeFeatureLayer to disable fields prior to merging.


texts: How to use Chrome element inspector to inspect marker icons in Google Maps API application?

texts:
I have some json data that are loaded with google maps API map.Data layer, I want to inspect the marker icon element (as shown as stars in the attached picture), what's the trick to use browser element inspector to do that?

Currently no matter how I try, the inspector just doesn't select the hovered stars (json markers from the Data layer).

Update

Data Layer marker icon can be inspected directly from dom tree.
The selector for these icons is 'img[src*="data:"]'. Through this way, additional css styles can be applied to the icons to improve the UX.


texts:
I'm really stumped with the problem. I have a polygon layer with no intersects or holes and I need to remove some 'islands', as defined by attributes of a certain size within the layer.

It's a hard thing to explain, so refer to the picture below. The different colors represent different attributes within my layer. What I need to do is remove the orange islands of certain size (less than 5000m²) that do not share  border with an Orange, Green or Neon polygon.

I cannot simply select all of the orange polygons under a certain size because that would remove some polygons that make up large areas of polygons.

I can't dissolve the orange than remove polygons under a certain size because I will lose the individual polygon shapes which are extremely important.

I am really stumped here guys, I'll take a simple solution or any advice. I am pretty handy with python so if you think a programming solution is best let me know.


texts: Overlay Union Geopandas improve performance

texts:
How to speed up the "union" of two shapefiles using geopandas?

I have two large shapefiles:

geodataframe_left.shape = (3610, 12) (500 MB) (link)

geodataframe_right.shape = (16396, 3) (200 MB) (link)

One file contains global province boundaries (GADM), The other shapefile contains Hydrobasin level 6. I've tried multiple approaches (PostGIS, Google Bigquery, ArcGIS) but nothing beats geopandas' one line command:

gdf_union = gpd.overlay(gdf_left, gdf_right, how='union')

there is however one problem. The command is extremely slow on larger datasets including mine. I've waited for almost an hour now without success. I'm running my code on a large amazon EC2 instance (32GB RAM).

How can I improve the performance of this operation? Is it possible to use multiprocessing for example or can I split the problem into smaller chunks so I can monitor progress?

UPDATE:

I tried updating to geopandas 0.4.0 .. still no result after running for two hours.
Used the use_sindex = True flag. Deprecated since 0.4.0 (always uses sindex)
I implemented the tile approach and found weird behavior. Details below.

I created a fishnet grid of polygons of 10 x 10 degrees that I use to clip both layers before performing the union. The processing times per cell differ wildly and there are 4 cells that take > 10 minutes to process. One cell takes even 40 minutes.

When the process runs in parallel I get the following error (script continues though):

AttributeError: 'IndexStreamHandle' object has no attribute '_ptr'
Exception ignored in: <bound method Handle.__del__ of <rtree.index.IndexStreamHandle object at 0x7f007cb42e10>>
Traceback (most recent call last):
  File "/opt/anaconda3/envs/python35/lib/python3.5/site-packages/rtree/index.py", line 875, in __del__
    self.destroy()
  File "/opt/anaconda3/envs/python35/lib/python3.5/site-packages/rtree/index.py", line 863, in destroy
    if self._ptr is not None:

What's interesting is that these cells are at the 180 meridian. I've uploaded the GPKG here

Maybe shapely doesn't like this hemisphere crossing stuff.

UPDATE 2:

Running ArcMap locally, the process takes appr. 10 minutes. Instead of using geopackages, I used shapefiles. file1 file2

Executing: Union "hybas_lev06_v1c_merged_fiona_V04 #;gadm36_1 #" C:\Users\Rutger.Hofste\Desktop\werkmap\union_benchmark\output\union_arcgis_global.shp ALL # GAPS
Start Time: Fri Nov 30 10:32:04 2018
Reading Features...
Processing Tiles...
Assembling Tile Features...
Succeeded at Fri Nov 30 10:43:02 2018 (Elapsed Time: 10 minutes 57 seconds)

I also tried "Union" (both default and SAGA) in QGIS but this process is also ridiculously slow. I ran for 30min and the progress bar was at 4%.

UPDATE 3:
I'm implementing the tiled approach however for a small number of polygons the result of an intersection in shapely is not a multipologon or polygon but a geometrycollection with multiple geometries including lineStrings or just LineStrings. I would expect the result of an intersection to be polygon or multipolygon.


texts: MapFish Print 3 - Multiple maps on a file

texts:
I'm using MapFish to print maps and i'm stuck on a problem : printing 2 maps on the same page of a PDF.

If i want to print one map i do that :

In Config.yaml :

[...]
attributes:
  map: !map
    maxDpi: 400
    width: 390
    height: 510
  title: !string {}
  [...]
processors:
- !reportBuilder # compile all reports in current directory
  directory: '.'
- !createScalebar {}
- !createMap
  inputMapper: {map: map}
  outputMapper: {mapSubReport: map}
  [...]

In my template i have something like :

[...]
<subreport>
    <reportElement x="0" y="40" width="390" height="510" uuid="fa145068-76a5-4834-98ed-ce65b1976b3d">
        <property name="local_mesure_unitwidth" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.width" value="px"/>
        <property name="local_mesure_unitheight" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.height" value="pixel"/>
    </reportElement>
    <subreportExpression><![CDATA[$P{map}]]></subreportExpression>
</subreport>
[...]

And the JSON i send :

"map": {
  "center": [
    148641,
    168747
  ],
  "scale": 1234,
  "layers": [...],
  "projection": "EPSG:31370",
  "dpi": "72"
}

That works perfectly.

So i tried to adapt to have 2 maps, one on the left one on the right :

In Config.yaml :

[...]
attributes:
  mapLeft: !map
    maxDpi: 400
    width: 390
    height: 510
  mapRight: !map
    maxDpi: 400
    width: 390
    height: 510
  title: !string {}
  [...]
processors:
- !reportBuilder # compile all reports in current directory
  directory: '.'
- !createScalebar {}
- !createMap
  inputMapper: {map: mapLeft}
  outputMapper: {mapSubReport: mapLeft}
- !createMap
  inputMapper: {map: mapRight}
  outputMapper: {mapSubReport: mapRight}
  [...]

In my template :

[...]
<subreport>
    <reportElement x="0" y="40" width="390" height="510" uuid="fa145068-76a5-4834-98ed-ce65b1976b3d">
        <property name="local_mesure_unitwidth" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.width" value="px"/>
        <property name="local_mesure_unitheight" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.height" value="pixel"/>
    </reportElement>
    <subreportExpression><![CDATA[$P{mapLeft}]]></subreportExpression>
</subreport>
<subreport>
    <reportElement x="400" y="40" width="390" height="510" uuid="fa145068-76a5-4834-98ed-ce65b1187d52">
        <property name="local_mesure_unitwidth" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.width" value="px"/>
        <property name="local_mesure_unitheight" value="pixel"/>
        <property name="com.jaspersoft.studio.unit.height" value="pixel"/>
    </reportElement>
    <subreportExpression><![CDATA[$P{mapRight}]]></subreportExpression>
</subreport>
[...]

And the JSON  :

"mapLeft": {
  "center": [
    148641,
    168747
  ],
  "scale": 1234,
  "layers": [...],
  "projection": "EPSG:31370",
  "dpi": "72"
},
"mapRight": {
  "center": [
    148641,
    168747
  ],
  "scale": 1234,
  "layers": [...],
  "projection": "EPSG:31370",
  "dpi": "72"
}

That doesn't work at all. Mapfish gives errors when loading the config.

So i guess i did it wrong, but i don't know how to make it work.

Does anyone have a suggestion ? Or even better a solution ?


texts: Having issues removing white background from Raster images

texts:
I've got an issue with 4 raster layers that I have converted from kml files into layer files. It would only allow me to convert them if I used ground overlay and now the 4 raster images are overlapping each other. I can't come up with a way to let them all be seen without the background obscuring the other images behind it. I've tried using symbology options and it doesn't do anything.

I've attached a screen capture. Basically the white background is what I'm trying to get rid of. I'm not sure if it's possible though.

Also I have the original PNG's and when I import the png's themselves and edit the background I can eliminate the white background however they have no spatial reference. When I try the same thing with these layer files you see in the screen cap, it doesn't work. Any thoughts?

Thanks!
-Brett


texts: QFileDialog remains open after selection made

texts:
I have created a QGIS plugin to import a tab separated value file, using:

filename = QFileDialog.getOpenFileName(self.dlg, "Select file to import ", path, '*.*')
uri="file:///"+filename+"?"+"delimiter=%s&xField=%s&yField=%s"%("\t", "Set Grid Easting", "Set Grid Northing")
vlayer = QgsVectorLayer(uri, "SpreadTests", "delimitedtext")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

It works as expected except that the dialog box to select the file stays open after a selection is made.  Is there something I am doing wrong? or something else I need to add?


texts: Using ModelBuilder to Iterate and extract subdatasets with in hdf file and project subdatasets?

texts:
I have a large file of n-files.  I am interested in using the model builder within GIS to extract subdataset fields within a Modis hdf file.  I am particularly interested in extracting Aerosol data from the hdf packed file.  I have to do this over several days (to years) worth of data.  Once I have extracted these subdatasets within the hdf file, I would need to then project the subdataset into WGS1984.

I am interested in the model building approach, so I can learn its functionality.  I know python may provide another solution, however, for the moment, I am interested in sticking with ModelBuilder (for educational purposes...).

I have attempted to create a model, but due to limited experience with this, I am confused as to why it is erroring out.  I thought as long as I iterate within all the hdf files within my directory, and pull the subdatasets, then project those subdatasets, I would have usable layers.  This is not the case.

I am using ArcGIS 10.2 for Desktop.


texts:
Is it possible to automatically refresh the QGIS map canvas in the background?
I have a .geojson file with the data. This data changes (in case of a simulation) but the map canvas only refreshes after I press the F5 button (refresh). In my case it should refresh the map every 100msec.
In the "layer properties/Rendering" there is a option "Refresh layer at interval (seconds)". I activated this but it has not the function I've searched for.
I found the command for the Python Console: iface.mapCanvas().refreshAllLayers().
This is exactly the function I've searched for. But I need it called periodically in the background.
I tried this with a self-written plugin, but without success. (Python script with an endless loop -> GUI freezes)


texts: Generating lines from coordinate file using R?

texts:
I have got an excel sheet with X and Y coordinates, actually Y is constant, since it is a straight line. I want a Shapefile with as many features as distance between points. I only get one line from the whole transect by using Linestring, but I need every segment, I tried with Multilinestring, but sincerely, I'm completely a beginner and after 2 days trying I couldn't get it.

library(sf)

x<- c(P$x)
y<- c(P$y)
coord2<- cbind(x,y)
print(coord2)
head(coord2)
L1<- st_linestring(coord2)
plot(L1, axes = TRUE)
class(L1)
geometria7 <- st_sfc(L1)
plot(geometria7, axes = TRUE)
st_write(geometria7, "proof.shp")


texts: Some characters in OpenStreetMap display as rectangles

texts:
I installed tileserver based on the instructions provided here for Ubuntu 14.04 using this approach . All fonts and Shapefiles were downloaded and installed from INSTALL.md.

The tileserver works well, but I have rectangles instead of some characters (in Chinese, Japanese, Korean, Indian, some African countries and cities names).

There is only debug information in log:

//initialization, no errors
//...
//loading fonts, no errors
//all fonts are loaded
DEBUG: Got incoming connection, fd 10, number 3
DEBUG: Got incoming request with protocol version 2
DEBUG: Got command RenderPrio fd(10) xml(default), z(7), x(113), y(51), mime(image/png), options()
//and so on.

This part of log of tile generating process (China tiles).


texts: Zoom Change Event Listener with google maps api

texts:
I'm trying to create an event listener that turns off a google fusion layer once the map is zoomed in to a certain level. I'm pretty new to javascript and its google maps api, so any help is appreciated. I found some code that seems close to what I need:

google.maps.event.addDomListener(map, 'zoom_changed', function() {
  var zoom = map.getZoom();
  if (zoom >= 1 && zoom <= 9) {
    layer_0.setMap(map);
    layer_1.setMap(map);
  }
  else if (zoom > 9 && zoom <= 11) {
    layer_0.setMap(map);
    layer_1.setMap(map);
  }
  else if (zoom > 11) {
    layer_0.setMap(null);
    layer_1.setMap(map);
  }
});


texts:
I've eaten at this amazing Brazilian BBQ restaurant in Calgary and I noticed its not on OpenStreetMap (OSM).

How do I add a label to OSM?

I am currently in the process of learning how to use QGIS and have been introduced to ArcGIS, I just need a refresher and some pointers on where to access reliable databases and upload my own layers or embed them into OSM [not sure exactly how its done at this point].


texts: How to get an ECW file into a database which QGIS can read?

texts:
I have ECW imagery I want to store in my spatial database. I am using QGIS and a PostGIS database. How do I do it?
Everything I have read so far (its actually unreadable to me) seems way too complicated when I was hoping for a drag'n'drop, cut'n'paste or "import to database" kind of solution. I'm a geographer not an IT boffin.
Anything out there, or should I buy an Esri license?


texts:
I'm configuring a server with several services and I want to know if is possible to change default Geoserver username/password using CLI, using configuration file/environment variables/CLI command/whatever, without having to log in manually and change it using the GUI.

Is this possible?


texts: GDAL .tif to MBtiles not creating all zoom levels

texts:
Using the following commands to convert a raster to an mbtiles file using gdal and then creating the ZXY schema. The zoom levels it creates are 2,3,4,5 but I need to be able to specify which zoom levels it creates, specifically zoom levels 9 through 12 for my leaflet map.
gdal_translate -of mbtiles temp.tif temp.mbtiles
gdaladdo -r nearest temp.mbtiles
mb-util --image_format=png temp.mbtiles ./temptiles

I've also tried setting zoom level like this, but then that does not produce any zoom levels, when extracted to the ZXY format via mbutil.
gdaladdo -r bilinear -oo "ZOOM_LEVEL=9" temp.mbtiles

Background info on the temp.tif original file:
Size is 1799, 1059
Origin = (-2699020.142521930392832,1588193.847443336388096)
Pixel Size = (3000.,-3000.)
Corner Coordinates:
Upper Left  (-2699020.143, 1588193.847) (134d 7'17.14"W, 47d50'44.64"N)
Lower Left  (-2699020.143,-1588806.153) (122d43'44.80"W, 21d 7'19.89"N)
Upper Right ( 2697979.857, 1588193.847) ( 60d53'28.48"W, 47d50'57.51"N)
Lower Right ( 2697979.857,-1588806.153) ( 72d16'48.48"W, 21d 7'28.62"N)
Center      (    -520.143,    -306.153) ( 97d30'21.52"W, 38d29'50.09"N)
Band 1 Block=1799x1 Type=Float64, ColorInterp=Gray
  Description = 2[m] HTGL="Specified height level above ground"
  Overviews: 900x530, 450x265, 225x133, 113x67, 57x34


texts: GDAL_Translate not compressing when NoData value is set

texts:
I am trying to compress a .TIF file using GDAL_Translate or GDAL_Warp. My Raster is 32mb to start with. When I use GDAL_Translate with the COMPRESS=PACKBITS option, it compresses to 1.9mb.

However if I now set the NoData value to -1 along with the COMPRESS=PACKBITS option my raster does not compress. Does anyone know what my be causing this? I have tried a number of - and + NoData values with the same result.