texts: FME Repeat features when a record matches an expression

texts:
I am working on automating the production of a land/parcel spatial table based on a state supplied cadastre (property layer) and injects from internal mssql systems. See this question for some further background FME Many to One to append attribute and not create new feature

When we have a parcel that is 0BUP48937 or 0SP117963 the SQL extract will have 1SP117963, 2SP117963 etc to match the number of different properties that exist within the parcel.

I currently match the state cadastre based on lotplan which then only incorporates the parcel with 0SP117963. How can I get FME to follow a logic like below

a - Where Lotplan starts with 0, match to the data from MSSQL (ex
   0SP117963)

b - then rather than rejecting the ones that start with 1,
   2, 3 etc (and have the same plan no) create a new polygon the same as 0SP117963

c -set the
   lotplan as 1SP117963 and then inject the rest of the relevant
   attribution.

So the table will be something like

The dataset does have Lot_num and Plan_num separately so the search for all records which have Lot_num=0 could also work...happy to also try a different logic if anyone can suggest something better.


texts: Renaming file geodatabase and keeping all mxd project linked?

texts:
I have a big folder (20GB) (with different folders with shapes, rasters and filegeodatabases)

This folder has a orthographic mistake that I would like to correct.

This gdb is involved in a big number of mxd projects so it is quite impossible to replace the path manually.

I have a way to move feature class one by one (and each of them take a lot of time) but this is too much work.

(I use arcgis 10.2)

Anyone has an idea to resolve my issue?


texts:
I can display a shapefile in openlayers but I need to add some attribute data.

I can open the dbf file in excel but there is no save as option in excel to the updated file as a dbf.

What is the best way (or software to use) to add attribute data to my shapefile?


texts: Is there a workaround for GRASS i.maxlik 'class' being a keyword in Python?

texts:
I am running i.maxlik in python - attempting to automate a processing chain.
I use the command:

grass.run_command('i.maxlik', group= group, subgroup=subgroup, class='classification', sigfile = sigfile, verbose = True)

Infuriatingly, however 'class' is a python keyword and it doesn't run due to 'invalid syntax'.

Would you know a way to work around this? I have tried emailing the GRASS developer list but no response.

Regards

Becky


texts: Expected traceback not shown on QgsRunProcess

texts:
On QGIS 2.18 python console (windows 10) I execute:

from qgis.core import QgsRunProcess
myRun = "C:\Python27\python.exe C:\Users\MyUser\Desktop\hw.py"
QgsRunProcess.create(myRun,True)

where hw.py is a simple python script that contains:

try:
    print("Hello world!")
    os.mkdir("aaa")
except Exception as e:
    print("My Exception has been caught: {} ".format(e))
os.mkdir("bbb")

As there's no import for os, the expected output in the dialog should be as follows (as it is when I execute in a windows console):

Hello world!
 My Exception has been caught: name 'os' is not defined
 Traceback (most recent call last):
 File "./Desktop/hw.py", line 6, in <module>
    os.mkdir("bbb")
 NameError: name 'os' is not defined

But what I get instead is:

Starting C:\Python27\python.exe C:\Users\MyUsers\Desktop\hw.py...
 Hello world!
 My Exception has been caugth: name 'os' is not defined

Done

with empty lines where the traceback should be shown.

Moreover I guess this could be related to stderrAvailable method (https://qgis.org/api/2.18/qgsrunprocess_8cpp_source.html  on line 103) but I can't see this error on the dialog, and I can't find how to get it shown there.


texts: Creating clusters of nearby points using QGIS?

texts:
I m trying to create clusters of nearby points (which would be the adresses of my customers around a store).
Buffer is not an option cauz i need the zone to be limited by customers, so a polygone shape not round
I've tried the convex hull but then my attibute table is empty and I need all the associated data inside the shape
I've also tried contour but i have 500 stores and i need to do it for each store separately
The point to polygone extension considers all the points together

Does someone have an idea of how to do this?

I'm working on QGIS 2.0.1.


texts: AutoDesk Infrastructure Studio - How to enable POST for WFS services?

texts:
I'm attempting to get feature info with the following url:

SERVICE=WFS&REQUEST=DescribeFeatureType&VERSION=1.1.0&TYPENAME=ns75382197:MAPSHEET_GRID

With many, many typenames listed after this first one.  Is there any way to enable the POST service analogous to this GET request?


texts: Why does QGIS 2.0 print composer not render the legend names?

texts:
I made a map using the print composer in QGIS 2.0. However, in the legend it does not display the name of all the items listed. I have updated the legend items, removed the items and added them again but does not work. Also, if I add an image to the map it is not displayed (item send to front does not make a difference).
Am I missing a setting or is it a bug?


texts: Using Spatial Analyst Classification in ArcGIS Engine?

texts:
I tend to use ML classification of Spatial Analyst ToolBox in ArcGIS Engine, I wrote this sample code:

        ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
        IAoInitialize aoInit = new AoInitialize();
        esriLicenseStatus licStat = aoInit.Initialize(esriLicenseProductCode.esriLicenseProductCodeEngine);
        licStat = aoInit.CheckOutExtension(esriLicenseExtensionCode.esriLicenseExtensionCodeSpatialAnalyst);

        Geoprocessor gp = new Geoprocessor();

        gp.AddToolbox(@"C:\Program Files (x86)\ArcGIS\Desktop10.1\ArcToolbox\Toolboxes\Spatial Analyst Tools.tbx");

        //--------------------------------------ENVI

        IVariantArray parameters = new VarArrayClass();
        parameters.Add(inputtxt.Text);
        parameters.Add(sigtxt.Text);

        parameters.Add(outputtxt.Text);

        IGeoProcessorResult result;
        result = (IGeoProcessorResult)gp.Execute("MLClassify ", parameters, null);

but unfortunately the code stop with this error:

Error HRESULT E_FAIL has been returned from a call to a COM component.

I do not know how to handle it?

The value of LicStat is         licStat esriLicenseCheckedOut

The Line      result = (IGeoProcessorResult)gp.Execute("MLClassify ", parameters, null);
is throwing the non-helpful error

As you said ---ArcGis 10 way is different and as that tool is a python script I'm not sure you can use it that way.--- what is its solution?, what should i do now?


texts: How to store complete edit history OSM PostgreSQL?

texts:
I successfully built a PostgreSQL server on an Ubuntu VM and populated it with the latest Geofabrik .osm.pbf complete dataset for Nepal.  Now, I'd like to append to that dataset all of the historic edit data.  I'd like to keep every version of every node, rel, and way, not just the latest, most up-to-date version.

My understanding is that OSM2PGSQL removes deprecated data when it appends the new data from a diff.  Is this correct?  If so, how can I go about storing a full history of edits, and not just the most recent ones?

When I try building an initial database with a full-history file like those found here, I get a crash with a duplicate key error, both on macOS and on Ubuntu.

  Processing: Node(2760k 32.1k/s) Way(0k 0.00k/s) Relation(0
  0.00/s)COPY_END for COPY planet_osm_ways FROM STDIN;  failed: ERROR:  duplicate key value violates unique constraint "planet_osm_ways_pkey"
  DETAIL:  Key (id)=(4748383) already exists. CONTEXT:  COPY
  planet_osm_ways, line 2

  Error occurred, cleaning up

I tried following the OSM-history tutorial found here, but the dependencies used are deprecated on the latest Ubuntu and 1) make it so you have to uninstall curl (!!) and 2) meant I would have to dig into some dependency stuff which is way about my pay grade in terms of my experience with Ubuntu.

To be clear, I just want to store the date stamp and author of every previous edit to any existing node, relation, or way.  I don't care what the actual edit was, I don't need to see it on a map.  I just want to know when it happened and who did it.


texts: Extracting multiple service polygons using ModelBuilder

texts:
I developed this model to perform multiple service area analyis. As you can see, in the final step I want to extract only the developed service polygons from the model and save them separatedly to work with. I do this by selecting only the polygons data and then copy these polygons to my disk.

For some reason I'm only ably to extract and copy the last developed service polygon, while I need all different service polygons for my further analysis. I tested my model and it works perfectly, untill this last step of extracting and copying the needed polygons. I do not know how I can adjust my model to extract and save each service polygon individually. In my previous models this problem was solved by using the percentage (%) symbols, but this doesn't seem to work here.

Is there an obvious thing I'm missing here?


texts: Making thematic layer 50% translucent in MapBasic?

texts:
I'm trying to make this enhancement to an old program written by another author.
I was able to make the layer translucent prior to Shade Window command, but its not working.
Here is the code:

Add Map Window nMapId Auto Layer TRADE_AREAS
Set Map Window nMapId Zoom Entire Layer TRADE_AREAS
Set Map Window nMapId Layer TRADE_AREAS
Set Map Window nMapId Layer TRADE_AREAS Alpha 127
    Label
        Position Center
        Font (Verdana, 0, 5, BLACK)
        With Zip_code
        Auto On Overlap Off

objBuffer = Buffer(objStore, 100, fRadius, "mi")
Insert Into MKT_BUFFER (obj) Values (objBuffer)

sCmd = "Shade Window " & nMapId & " TRADE_AREAS With Trade_Area Values "
        + """Primary"" Pen (1,2,0) Brush (2,16711680,16777215), "  'RED, WHITE
        + """Secondary"" Pen (1,2,0) Brush (2,65280,16777215), "   'GREEN, WHITE
        + """Tertiary"" Pen (1,2,0) Brush (2,16776960,16777215) "    'YELLOW WHITE
Run Command sCmd


texts: Coloring lines with polygon grid using ArcGIS for Desktop?

texts:
In Arcgis, I have a polygon grid layer (100 cells) and a layer with polylines (the polylines cross the polygon grid). The polygon grid layer has in the attribute table a numeric integer field (named "weight") that give to the respective cell a weight. Example:
cell1-->100
cell2-->105
cell3-->2
cell4-->27
cell5-->3
...-->...
cell100-->89

my goal is ONLY purely visual nothing more. I want color each cell of the the polygon grid with a color based on the weight attribute (clearly if weight of cells is almost equal I want use similar color). After (this is the main problem) I want that only a part of each polyline take the color of the cell that cross. So the final view of each polyline will be N-polyline with multiple gradient color. The intersect tool isn't a good solution because in real case I have a lot of polyline and a lot of cells, so the intersection process is very slow.


texts: Calculating area of newly added polygons on same layer?

texts:
I have a layer where I had 6 polygons. I then used the field calculator so I could get a column with the area of each of these polygons and I thought this would be applied to any other polygons I made after that, which was not the case since every new polygon appears as "NULL" in the area column.
How do I make so the area is calculated automatically for all the rows?


texts: Hiding points in shapefile layer using QGIS

texts:
I am using QGIS 3.4

I would like to hide all points that end with "5". For example, I would like to hide: 55, 65, 75, 85 and any point ending in "5".

What method can be use to make this possible?


texts:
I'd like to make a map where you can click a state and have some information in the comment box, but I'd also like to include US territories (Guam, PR, USVI, etc.). These are pretty far from the continental US, so it's not super convenient or attractive to set the main map extent to include them.

Is there a tool in ArcGIS online (story maps, story map templates, web apps, etc) that would allow me to easily make buttons that can be clicked to provide comment box-like information?
Something like this: https://www.cdc.gov/coronavirus/2019-ncov/cases-updates/cases-in-us.html is what I have in mind.


texts:
I have the park boundary for Point Reyes National Seashore (https://irma.nps.gov/DataStore/Reference/Profile/2224545?lnv=True) and the NHD data for that park boundary. I want to do a drainage density analysis, but I don't want to include the 400m of water that extends beyond the shoreline in my area, not to mention the bay. I want to clip the park boundary to the extent of the shoreline, but the shoreline is an incomplete polygon (the park is connected to land!) so I can't actually clip it.

What is the best method to clip this to the seashore in a relatively quick manner?


texts: WFS 2.0 with get request outputformat gml3 returns body using wfs 1.1.0 schema, is this expected?

texts:
I am making GET requests against a wfs service specifying version 2.0, when I do not specify an outputformat in the request it returns a body using the schema namespace for wfs 2.0 (wfs:), when I specify outputformat=GML3 it returns a body using schema namespace for wfs 1.1.0 (gml:)

I know that in the get capabilities it states that it supports gml3.2, gml3, gml2 and that these are largely associated with wfs 2, 1.1 and 1 respectively but I was wondering whether the use of a previous wfs version schema in return body was expected behaviour in this scenario as I wasn't expecting it and I can't find any real info on this.


texts: Why do I get "Invalid Data Source" when opening a GPX file from ETrex 30?

texts:
Why is it if i load the .gpx files in qgis it says "Invalid Data Source" as well as in google maps?

I want to view the waypoints from the .gpx file. I used gpx etrex 30


texts: Survey123 v 3.12.274 no longer submits in background

texts:
I have a collector-survey123 integrated workflow on iOS for field personnel. They select a route and launch a survey from a custom url in the collector pop-up. This pre-populates some questions etc, and when they submit the survey it brings them back to collector with a callback.
Since the last update, when the callback fires they are brought back to collector, but it seems that the surveys are no longer sending in the background, and so just sit in the outbox or are not sent until they launch the next survey instance, at which point the submit finishes and then they are brought to an empty outbox, meaning they need to return to collector and re-launch the survey.
It does not happen every time, but most now. It seems like the survey is not fully sending before the callback fires, and survey123 is no longer submitting/running in the background. Any suggested workaround? I changed the callback to callback:submit, but that seemed to have no impact.


texts: Merging Sentinel 2 RGB bands with rasterio

texts:
I'm trying to merge Sentinel 2 bands RGB bands with rasterio like this:

import rasterio
import numpy as np
from rasterio import plot
from rasterio.plot import show

band2=rasterio.open("B02.jp2")
band3=rasterio.open("B03.jp2")
band4=rasterio.open("B04.jp2")

rgb=rasterio.open('rgb.tiff', 'w', driver='Gtiff',
                          width=band2.width, height=band2.height,
                          count=3,
                          crs=band2.crs,
                          transform=band2.transform,
                          dtype='uint16')
rgb.write(band4.read(1),1)
rgb.write(band3.read(1),2)
rgb.write(band2.read(1),3)
rgb.close()

This works as expected and creates a RGB-GeoTIFF that is properly displayed in QGIS. However if I try to visualize it with rasterio, I get an almost black image:

show(plot.adjust_band(rasterio.open('rgb.tiff').read([1,2,3])))

when I inspect the normalized values the mean for each band is around 0.05. Is there something wrong how I create the RGB?


texts: Why do two tables with valid geometry return self-intersection errors when intersected?

texts:
I have two tables A and B, both have valid geometries when I check with ST_IsValid. In my query, I intersect the two and the query stops in the middle because of this error:
ERROR:  lwgeom_intersection: GEOS Error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point

I read ST_Union fails with TopologyException despite valid polygons and using ST_SnapToGrid which seems very similar, and tried the ST_CollectionExtract as mentioned in the accepted answer, but I still the same error.

When I run the query on a subset of the data I notice some records are ST_MultiPolygon while others are ST_Polygon, is this causing the issue and, if so, how can I fix it?

The original query has other fields which I've left out here but it generally looks like this:

SELECT
   ST_CollectionExtract(ST_Intersection(ST_Transform(r.geom, 5070),
                                        ST_Transform(c.geom, 5070)),3) as geom,
   a.place,
   b.name

FROM table1 as a
JOIN table2 as b
WHERE a.geom && b.geom
AND ST_Intersects(ST_Transform(r.geom, 5070), ST_Transform(c.geom, 5070))

I am projecting it because I also calculate areas in the full query.


texts: Tracing vector map from GeoTIFF for Garmin

texts:
My goal is to create a vector map for my Garmin GPS by tracing over a topo raster.  The vector data would be exported in ESRI SHP format and later converted for Garmin.  So far, I can load the geotiff layer, however am stumped on creating a vector layer.

Where do all the shapes come from to trace on the map?


texts: Packaging succeeded, but publishing failed to ArcGIS for Server?

texts:
I am using ArcGIS Server 10.1. I have tried all the suggested solutions above this thread but still unable/failed to create/execute service. The error message I get is:

  Packaging succeeded, but publishing failed. ERROR 001369: Failed to
  create the service Please consult the GP Results Window for enhanced
  error description


texts: "Runtime error: cannot open workspace" after creating a connection file and changing the version

texts:
I'm trying to run an update cursor. Usually if I have an issue with databases with enterprise connections, making a connection file clears it up. This time when I tried to execute my code I just get a Runtime Error: cannot open workspace message.

I went ahead and changed the version but that didn't work either. I also tried using the workspace without the connection file and I got the same error. I've yet to make it past edit = arcpy.da.Editor(workspace)

Here's the bulk of my immaculate code:

def build_connection_file():
    """This function handles everything involved in building a connection
file to the feature class"""

    directory = r"C:\TempConnectionFiles"
    if not os.path.isdir(directory):
        os.path.makedirs(directory)

    feature_dataset = r"ThisIsReal"
    feature_class = r"ThisIsReal"

    sde = r"Database Connections\Real.sde"

    versions = [version.name for version in arcpy.da.ListVersions(sde)\
                if (getpass.getuser() in version.name)]

    print versions[-1]

    connection_name = "ThisIsReal.sde"
    if not os.path.isfile(os.path.join(directory, connection_name)):
        arcpy.CreateDatabaseConnection_management(
            out_folder_path = directory,
            out_name = connection_name,
            database_platform = "real_platform",
            instance = "real_instance",
            account_authentication = "OPERATING_SYSTEM_AUTH",
            username = "",
            password = "",
            save_user_pass = "DO_NOT_SAVE_USERNAME",
            database = "real_database",
            version_type = "TRANSACTIONAL",
            version = versions[-1])

    connection_feature_class_full = os.path.join(
        directory, connection_name, feature_dataset,
        feature_class)

    print connection_feature_class_full

    #this was required last time I did this.
    table_view = arcpy.MakeTableView_management(
        in_table = connection_feature_class_full,
        out_view = r"connection_table", where_clause = "",
        workspace ="", field_info = "")

    arcpy.ChangeVersion_management(
        in_features = table_view, version_type = "TRANSACTIONAL",
        version_name = versions[-1], date = '')

    print("changed version")

    return connection_feature_class_full

connection_feature_class = build_connection_file()
edit = arcpy.da.Editor(connection_feature_class_full)
edit.startEditing(True, False)
edit.startOperation()

Is there a setting in my database that could be responsible for this? Is there something wrong with my code?


texts: Setting Map Algebra/Raster Calculator's pixel depth

texts:
I use Raster Calculator a lot for my work and was wondering if there is a way to set the output pixel depth to a desired value without using Copy Raster tool after I run a calculation, which is the only facility to do that.

In the ArcGIS 10.x help page, it mentions about the mechanism which controls how the no data value is created and promotion of pixel depth to the next higher bit, if required. However my experience with raster calculator is that whenever I use float rasters, it yields a 32-bit raster, even though the value range is limited to 8-bit signed or even binary. The best example to this is NDVI, even though Calculate Statistics tool says it ranges between float -1 and +1, I have this saved as 32-bit.


texts:
I'm trying to create a buffer on a projected file and the result is a huge blob.
I've read that this can be linked to the CRS used, so I made sure that I'm using a projected CRS (EPSG:32734 - WGS 84 / UTM zone 34S - Projected). The buffer specifies meters (not degrees like when the CRS is not a projection).
I'm running out of ideas ...
Windows10 64bits, qgis 3.14.16-Pi
Roads.gpkg


texts: Moving gdal2tiles.py output directly to s3 from inside a fargate task

texts:
I am running the latest version of gdal in a docker container built using:
FROM osgeo/gdal. The docker container is running on AWS as a fargate task.
I want to use gdal2tiles.py to create a large tileset that covers all of Germany.
The tileset is generated from 52 sentinel images in the jp2000 format. First I warp all files that are not in UTM32 (EPSG:32632) into UTM32 using:
gdalwarp -r cubic -t_srs EPSG:32632 {file} {warped_file}
After I have all the files I need in the same projection, I build a vrt containing all 52 files:
gdalbuildvrt -o { vrt_name} {os.path.join(directory, '*.jp2')} -srcnodata 0 -te {xmin} {ymin} {xmax} {ymax}
After that I use gdal2tiles.py like this:
gdal2tiles.py --zoom=17 --processes=4 --webviewer=none { vrt_name} {tiles_dir}
The tiles are created correctly, my problem is that the fargate task has only 20 gb of ephemeral storage. That storage is not enough to store the thousands of PNGs created by gdal2tiles.py. I know the only other way to increase storage for a fargate task is to mount an EFS volume to the task. However, for budgetting reasons and because the tiles will get served from an s3 bucket, I do not want to use EFS. Is there a way to move the output of gdal2tiles.py directly to s3, without storing it locally?
I have read about virtual filesystems for gdal here:
https://gdal.org/user/virtual_file_systems.html#vsis3-aws-s3-files
But this seems to only work for reading files or writing sequential files to s3. Or am i mistaken there?
Does anyone know a solution or workaround?


texts: Use OpenLayers 6 node package with jQuery

texts: I want to use ol node package from nodejs. This package is written in ES6. I am using Asp.net 4.5(aspx) and jQuery to use this. I would like to know how I can support multiple tile layer source with ASPX and jQuery, is there pseudocode available.


texts: Adding custom map tile layer using QGIS?

texts:
I am trying to display the following background map in QGIS 2.0:
http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/${z}/${y}/${x}.
In a custom web map client, this is quite easy:
http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/xyz-esri.html

I tried to use the "OpenLayers plugin" in QGIS, but there the list of background tile maps is predefined and I cannot add my own custom map tile layer.

Is there any easy workaround (without any programming) how I could show the custom map tiles layer in QGIS?


texts:
I am beginner in Linux. I have a problem to use Shapely 1.6.4 post2...
I installed shapely by this command:

$ pip install shapely

and it worked good

Collecting shapely
  Using cached https://files.pythonhosted.org/packages/81/d1/b8e1b089a8ddd6df74be583d70373eac55c725c6197c115efbd3c3e1509f/Shapely-1.6.4.post2-cp27-cp27mu-manylinux1_x86_64.whl
Installing collected packages: shapely
Successfully installed shapely-1.6.4.post2

but when i want to import shapely...

$ python chp1.py
Traceback (most recent call last):
File "chp1.py", line 3, in <module>
from shapely.geometry import Point
ModuleNotFoundError: No module named 'shapely'

I use python 3.6.7

Could you explain the problem?


texts: How can I get a join for pre-existing fields?

texts:
I am joining some values from a table to a shapefile. The table and shapefile have the same column names and types, but when I join, Arcmap creates a new field, and I just want the join to write the data from the table into the fields that already exist.

How would I do that?


texts: Getting the index of a node in a qgslayertree

texts:
Can I get the index ([integer]) - not the ID- of a child in a QgsLayerTreeGroup ?

I'd like to insert a layer using insertChildNode(int index, QgsLayerTreeNode *node) in the layer tree. index is required for control of the exact position of the inserted layer in the LayerTreeGroup.
The exact position of the inserted layer in the layer-tree depends on other layers. So i can not give a default value as index.  I've to calculate this value depending on the index of the otherLayer. For this reason I'd like to determine the index-value of the otherLayer. In the API-doc I only read about functions to determine its ID - which seems not to be the same as the index.


texts: Converting every pixel of TIFF file from decimal to binary

texts:
I have a folder of .tif-files. I would like to convert every value in each of these files from decimal to binary. This code mostly works. The inner loop correctly outputs the binary equivalent for every value in the .tif-file. The outfile is also written to the output folder with the correct name and extension (.tif).
However, while something is being written to the output file (each one has a size of about 300 mb, comparable in size to the original .tif), the file cannot be opened in ArcGIS. When I try I get the error:

Failed to add data, unsupported data type.

How can I properly save the output in .tif format to be accessible in ArcGIS?
I'm using Python 3.6.1 in the Jupypter notebook in ArcGIS Pro 2.5.
import arcpy, os, glob
from arcpy.sa import *

arcpy.CheckOutExtension("Spatial")
arcpy.env.addOutputsToMap = 0
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.overwriteOutput = True

input_ws = glob.glob(r'D:\MODIS_DOWNLOADS\Australia_500\QA\*.tif')
output_ws = r'D:\MODIS_DOWNLOADS\Australia_500\QA\BINARY'

for r in input_ws:
    basename = os.path.basename(r)
    r1 = arcpy.Raster(r)
    NUM_COLUMNS = arcpy.management.GetRasterProperties(in_raster=r1, property_type="COLUMNCOUNT", band_index="")[0]
    COLUMNS_1=int(NUM_COLUMNS)
    print(COLUMNS_1)
    NUM_ROWS= arcpy.management.GetRasterProperties(in_raster=r1, property_type="ROWCOUNT", band_index="")[0]
    ROWS_1=int(NUM_ROWS)
    print(ROWS_1)
    outname = "BINARY_"+basename
    print(outname)

    outfile= open(os.path.join(output_ws, outname),"w+")

    try:
        for i in range (COLUMNS_1-1):
            for j in range (ROWS_1-1):
                if r1[i,j]>0:
                    H=(r1[i,j])
                    BINARY_RAW=bin(H)
                    BINARY_NEW=BINARY_RAW[2:]
                    outfile.write(BINARY_NEW)
    finally:
        outfile.close()

Similarly, I've tried solving this problem using the raster cell iterator, available beginning in ArcGIS Pro 2.5. This code runs and throws no errors, saving an output .tif file to the output folder. However, while the output file opens in ArcGIS Pro, the file is entirely made up of NaN.
import arcpy, os, glob, numpy
from arcpy.sa import Raster, RasterCellIterator

arcpy.CheckOutExtension("Spatial")
arcpy.env.addOutputsToMap = 0
arcpy.env.parallelProcessingFactor = "100%"
arcpy.env.overwriteOutput = True

ws1 = glob.glob(r'D:\MODIS_DOWNLOADS\Australia_500\QA\*.tif')
outws = r'D:\MODIS_DOWNLOADS\Australia_500\QA\BINARY'

for r in ws1:
    basename = os.path.basename(r)
    r1 = arcpy.Raster(r)
    outname = "BINARY_"+basename
    print(outname)
    outRas=Raster(r1.getRasterInfo())
    outRas.readOnly=False

    with RasterCellIterator({'rasters':[r1,outRas]}) as rci:
        for i,j in rci:
                if r1[i,j]>0:
                    H=(r1[i,j])
                    BINARY_RAW=bin(H)
                    BINARY_1=BINARY_RAW[2:]
                    BINARY_2=int(BINARY_1)
                    outRas[i,j]=BINARY_2
        outRas.save(os.path.join(outws, outname))


texts: Semantics of polygon points - How to join the dots?

texts:
The question arise from reading the standard "OGC® Sensor Planning Service Interface Standard 2.0 Earth Observation Satellite Tasking Extension", that you can get from http://www.opengeospatial.org/standards/sps (Although, I think it is more general and also applies to kml standard too).

Suppose I receive a Polygon expressed in terms of WGS84 UTM latitude longitude points, how should I interpret that set of points? I mean, how should I build the polygon from the points?:

Should I use straight lines in the lat lon plane (linear interpolation between points in the lat lon plane)?
Should I use great circle curves to join contiguous points?
Should I use the geodesic curves in the WGS84 ellipsoid?
Must this be explicitly specified in the interface?

Specifically, I'm asking about the RegionOfInterest class that contains a Polygon class in page 31 of that standard.


texts: How to create rectangles of varying dimensions from a multi-point Feature Collection in Google Earth Engine?

texts:
I have a Shapefile containing point data that has three variables (i.e. name, width, height). I ingested it as a Feature Collection in Google Earth Engine (GEE). I want to make rectangles of varying dimensions using the 'width' and 'height' data in GEE. I do not know how to start doing this. In GEE, the closest thing I can do is to create square buffers around a point using .buffer() and .bounds() as I show in the sample code below but that is not what I want.
Can anyone advise on how I can code to create geometries of varying dimensions from point data in GEE?
I have made this GEE asset publicly available for anyone to use.
This is the attribute table of the ee.FeatureCollection.

Here is the sample code I mentioned that creates a fixed distance square buffer around each point. I am not sure if this can be modified to create rectangles of varying widths and heights according to the values in the attribute table.
//  Load point data
var points = ee.FeatureCollection('users/earthorbgis/GEE_help/gee_point_sample');

//  Buffer points to 500 x 500 m squares
var buffer_interval01_MMR = points.map(function(f) {
  return f.buffer(250).bounds();
});

//  Create custom rectangles using point data columns Width and Height
//  Unsure how to start ...


texts: QGIS Spatialite error: No Such Function: spatialite_version

texts:
I am trying to use Midvatten on QGIS 2.18.13 on a mac.  When I create a new database, I get

  "no such function: spatialite_version" (see error message below).

I've uninstalled and reinstalled QGIS (using Kyngchaos's packages) and Midvatten. None of that has helped.

How do I get this basic spatialite function working?

  Traceback (most recent call last):

  File "/Users/Becca/.qgis2/python/plugins/midvatten/midvatten.py",
  line 853, in new_db
      newdbinstance = newdb(verno)

  File
  "/Users/Becca/.qgis2/python/plugins/midvatten/tools/create_db.py",
  line 41, in init
      self.create_new_db(verno,user_select_CRS,EPSG_code, delete_srids)  #CreateNewDB(verno)

  File
  "/Users/Becca/.qgis2/python/plugins/midvatten/tools/create_db.py",
  line 82, in create_new_db
      versionstext = self.cur.execute('select spatialite_version()').fetchall()

  OperationalError: no such function: spatialite_version


texts: Convert CSV file with long and lat geometry column to sf object in R


texts: Upload many layers from PostGIS database to GeoServer using curl

texts:
I have more than 200 layers in my PostGIS database and I want to publish those layers in GeoServer using curl. I've imported one layer using the command below but can't figure out how to upload all layers in one step.

curl -v -u admin:geoserver -XPOST -H "Content-type: text/xml" -d "<featureType><name>testlayer</name></featureType>" http://localhost:8088/geoserver/rest/workspaces/myworkspace/datastores/mystore/featuretypes


texts: I create lines in JOSM Editor and while saving it to osm format it works well. saving and reopening it , but for saving it in osm.pbf format it works but in  reopening it, I get this error


texts: ArcGIS for Server: Get map service cache scales programmatically

texts:
I am trying to write a script to build tile caches for a bunch of map services. My issues is that to run arcpy.ManageMapServerCacheTiles_server() I need to pass a list explicitly stating scales at which to build tiles. Easy enough when the service is using the 'ArcGIS Online / Bing Maps / Google Maps' predefined tiling scheme because I know the exact scale values. But I would like to write this script so that it can be used for any map service, even if those with a custom tiling scheme.
To do this, however, I would need be able to lookup what scale values exist in the services tiling scheme. Is there any way to do this programmatically?


texts: Converting LiDAR point cloud to 3D mesh?

texts:
I have an unstructured point cloud of a street in LAS format (.las), scanned from a driving car. The point cloud also stores RGB information.

Is there any software that can convert the point cloud to a 3D mesh, for example FBX or OBJ? How?

I've tried:

Autodesk Recap Pro - it doesn't work with unstructured point clouds.
RealityCapture - it's not compatible with LAS format.
Pix4D - it requires photos for reconstruction.


texts: Setting 'classificationMode = Continuous' to a raster layer using PyQGIS?

texts:
I need to add to this script:

How to update raster color ramp with new min/max values using PyQGIS?

this function: QgsColorRampShader for setting 'classificationMode = Continuous', similarly as setColorRampType(QgsColorRampShader.INTERPOLATED), but I tried it and it doesn't work.

    def setStyle(self, layer):
    provider = layer.dataProvider()
    extent = layer.extent()

    stats = provider.bandStatistics(1, QgsRasterBandStats.All, extent, 0)

    if stats.minimumValue < 0:
        min = 0
    else:
        min = stats.minimumValue

    max = stats.maximumValue
    range = max - min
    add = range // 2
    interval = min + add

    colDic = {'red': '#ff0000', 'yellow': '#ffff00', 'blue': '#0000ff'}

    valueList = [min, interval, max]

    lst = [QgsColorRampShader.ColorRampItem(valueList[0], QColor(colDic['red'])),
           QgsColorRampShader.ColorRampItem(valueList[1], QColor(colDic['yellow'])),
           QgsColorRampShader.ColorRampItem(valueList[2], QColor(colDic['blue']))]

    myRasterShader = QgsRasterShader()
    myColorRamp = QgsColorRampShader()

    myColorRamp.setColorRampItemList(lst)
    myColorRamp.setColorRampType(QgsColorRampShader.INTERPOLATED)
    # TODO: add classificationMode=Continuos
    myRasterShader.setRasterShaderFunction(myColorRamp)

    myPseudoRenderer = QgsSingleBandPseudoColorRenderer(layer.dataProvider(),
                                                        layer.type(),
                                                        myRasterShader)

    layer.setRenderer(myPseudoRenderer)

    layer.triggerRepaint()


texts: What does "Place of Work only" mean in terms of US census ACS geographical restrictions on a dataset?

texts:
I'm trying to get data on commute methods for an area using American Community Survey 2015 table B08505A, but as you can see in this table describing geographic restrictions for ACS data: https://www2.census.gov/programs-surveys/acs/summary_file/2015/documentation/tech_docs/ACS_2015_SF_5YR_Appendices.xls
The data is reported for "Place of Work Only".  What does "Place of Work Only" mean in terms of US census geographic restrictions?  I can't find info on this anywhere.


texts:
I want to select an edge of parcel instead of creating multiple displacement links.
The problem is when selecting the parcel, and the edge of the parcel contains vertices, the links that will be attached to the parcel edge will cover the area between two vertices,
How to attache the links among all the parcel edge, regardless the number of the vertices.


texts: Convenient way to download MCD19A2 product using R

texts:
I want to download MODIS MCDA219 product. I have been downloading manually from

https://ladsweb.modaps.eosdis.nasa.gov/
and https://earthdata.nasa.gov .

As a new to R user, recently I have learn downloading this product using R with help of Download MODIS/MCD19A2 AOD product with R .
Problem is, the daily product are creating too much folder. For example if I try to download January month of 2018, then it will create folder for each day. So for January month, 31 hdf format files are in 31 folders each. How can I download MCD19A2 product by month and the HDF file?
I will need more than 50,000 dataset and I was looking for more convenient way to download the product using R


texts: Cross tab of two rasters across different polygons in a shapefile

texts:
I need to find the number of pixels of raster 1 that fall within different pixel values of raster 2. For example, there are 250 pixels of value=2 in raster1 that in value=89 of raster 2. Hence, I am using crosstab() from the raster package.

Now I need to find the cross tab of pixel count from the two rasters, within polygons in a shapefile. How do I do that?

One can think of this problem as- area (pixel count* pixel size) of opportunity (0-no opportunity, 1-opportunity specified in pixels in raster 1) in different categories of standards of living (1- low standard of living, 2-medium, 3-high specified in pixels in raster 2), across different US counties/states (polygons in shapefile).


texts:
I have a PostgreSQL/PostGIS table "geometry_table" that has a 1-to-many relation to another table "properties_table", where each geometry has many properties. These properties hold style values of geoemtries (e.g color, stroke-width, ... etc).
I have a CSS style document that's used to style this geometry_table layer, something like the following:

/*Lines*/
[GetStyleProperty(geomId, 'visible') = 'true'] [dimension(geometry) = 1] [@sd < GetStyleProperty(geomId, 'feature_max_scale')] [@sd > GetStyleProperty(geomId, 'feature_min_scale')]{
   stroke: [GetStyleProperty(geomId, 'stroke_color')];
   stroke-width: [GetStyleProperty(geomId, 'stroke_width')];
   stroke-dasharray: [GetStyleProperty(geomId, 'stroke_dasharray')];
   stroke-opacity: [GetStyleProperty(geomId, 'opacity')];
}

The function "GetStyleProperty" is a custom Geotool function I've created to read a given property name of a given geometry from the DB. This function first opens a DB connection, then executes a DB stored function, then closes the connection.

My problem is that loading the geometry_table layer takes too much time (more than 6 seconds to display arround 10 featrues only). I have recognized that the problem is from my CSS Style document, because if I use another style document (e.g the generic one) the layer is dispalyed in less than one second. I've measured the execution time of "GetStyleProperty" function and found that quarter of the time is wasted in opening the DB connection.

So, my question is: instead of opening a DB connection in every function call, is there any shared connection that I can use to connect to my DB? Or, how can I make a shared DB connection that lives as long as the GeoServer application lives?

I'm using the following code to create a DB connection:

public static Connection createConnection() throws ClassNotFoundException, SQLException {
        Class.forName("org.postgresql.Driver");
        String url = "jdbc:postgresql://localhost:5432/spatialdbnew";
        Connection conn = DriverManager.getConnection(url, "postgres", "postgres"); //use postgres username to skip row-level security to enhance performance

        /*
         * Add the geometry types to the connection.
         */
        ((org.postgresql.PGConnection) conn).addDataType("geometry",
                (Class<? extends PGobject>) Class.forName("org.postgis.PGgeometry"));
        ((org.postgresql.PGConnection) conn).addDataType("box3d",
                (Class<? extends PGobject>) Class.forName("org.postgis.PGbox3d"));

        return conn;
    }


texts:
I try to import a dwg file (using Projekt > DWG/DXF-Import) into QGIS (2.18.3 32Bit on Windows) which works quit nice.

Since I try to use the fid, I found that the numbers are not the same in QGIS as the one in AutoCAD MAP.

Is there a way to get the proper numbers directly in QGIS?


texts: How to accelerate v.generalize on massive vector data?

texts:
I have a big problem with my vector dataset (more than 1 million entities) and v.generalize in grass 7.0.5 and 6.4 (64 bits for both). I try to use Douglas method on this dataset but it is very very slow and I'm searching how to improve this treatment.

Firstly, i can't split my dataset for this job and I must conserve the topology.

For test, I use the lsat7_2000_80 raster from North Carolina data set for grass 7. Firstly, i vectorise (r.to.vect) this raster in areas with -s flag (smooth corners). I got 1642301 vertices for 224752 areas. Secondly, I use v.generalize with douglas method and treshold 15 meters. This step take 48 minutes (i5 2.6ghz, 8go ram on windows 7 64bits and Ubuntu 16.04 64 bits) and i search how accelerate it.

Does anyone have an idea on how to increase the speed of this tool?


texts: Locating the centroid (center of mass) of spherical polygons

texts:
I'm trying to work out how best to locate the centroid of an arbitrary shape draped over a unit sphere, with the input being ordered (clockwise or anti-cw) vertices for the shape boundary. The density of vertices is irregular along the boundary, so the arc-lengths between them are not generally equal. Because the shapes may be very large (half a hemisphere) it is generally not possible to simply project the vertices to a plane and use planar methods, as detailed on Wikipedia. A slightly better approach involves the use of planar geometry manipulated in spherical coordinates, but again, with large polygons this method fails, as nicely illustrated here. On that same page, @Cffk highlighted a paper which describes a method for calculating the centroid of spherical triangles. I've tried to implement this method, but haven't yet succeeded in calculating centroids that seem to fall in the proper place.

My attempt is implemented in Python and posted on StackOverflow here. I can copy the code over if wanted, but basically what I want to know is if the method described in the paper is applicable to a more complicated polygon boundary. Specifically by (1) stepping the base of the triangles along paired vertices (c.w. or c.c.w. ordered) of the polygon boundary, (2) calculating sub-centroids for each of these triangles and then (3) adding the area-scaled elements of the sub-centroids together to get a centroid value for the polygon as a whole.

In short, I am treating the equations in the paper as a 3d extension of the commonly-used method of determining the area and centroid of a planar polygon, but something doesn't appear to be holding up, and at the moment it seems more likely a flaw in the math or my logic than in the construction of the python code.


texts: Use American Polyconic Projection North American Datum 1927 with GDAL and QGIS

texts:
I am trying to use modern USGS DEM data with older USGS maps. The first step, as far as I am aware, is to reproject the contemporary DEM into the same projection as the old USGS map. The old USGS map notes that it's a "polyconic projection. 1927 North American Datum. 10,000-foot grid based on Oregon coordinate system, north zone."
Based off of some sources I found online, I need to define the projection using a proj code. A little more online digging yielded this:  "+proj=poly +lat_0=0 +lon_0=-122.625 +x_0=0 +y_0=0 +a=6377397.16 +b=6356078.96 +units=m +no_defs". I tried using this with gdalwarp, as follows:  gdalwarp -t_srs "+proj=poly +lat_0=0 +lon_0=-122.625 +x_0=0 +y_0=0 +a=6377397.16 +b=6356078.96 +units=m +no_defs" -r bilinear -of GTiff -co "TFW=YES" dem-edit.tif dem-edit-2-5.tif
The final image generated appears rotated, both if I open the TIF directly, and if I open it in QGIS.

My intuition is that this is something to do with the projection I defined.
Does anyone know how to correctly use American Polyconic, NAD 1927 with GDAL and CLI also QGIS? It seems that theres no EPSG code, and my best guess is that there's something wrong with the proj code.
Ultimately I plan to load the DEM and the old map into Blender for some 3D renderings, so if there is another projection that would give me close enough results to the original then that could probably work as well. In the GDAL command instead of defining a and b, I also tried +ellps =clrk66 to no avail.
If anyone has any advice or ideas, I'd be really grateful.


texts: geoserver/openlayers: multiple layer geotiff

texts:
I have a geotiff file with multiple layers (46 layers), each layer corresponding to a different date (the geotiff layers corresponding to sea surface temperature, 1 layer per year).

I imported this geotiff in geoserver, but I can't render it with the geoserver visualisation tool (based on openlayers).

1./ does geoserver support multiple layers > 3 bands?
2./ if yes, when I query the WMS in openlayers, how can I select a given layer?


texts: I am using LiDAR data to locate Sinkholes(depressions) in an area. I am doing this by making the contours and checking the contours(if they are circular enough or not, elongatedness, etc.). But I have a problem, I need to check the jaggedness of the circle edge, so If it is too jagged, I eliminate that contour. Any idea how can I do that?


texts: Choosing projection for calculating 2D 3D area in New Zealand?

texts:
I'm trying to use Surface Volume tool in ArcMap to calculate 2D and 3D area of a catchment.

What would be the best projection to use and how do i check the maximum error of the estimation?


texts:
I've downloaded us region boundaries from US census bureau, which is a .zip file. I can open to view the file with QGis as follows (4 regions + Alaska):

Now I have a dataset containing a lot of data point. Each data point has a pair of latitude and longitude. I want to use these numbers to identify which region a data point belongs to.

How can I do it with QGIS?


texts: Create points at specified distance from each other in random fishnet's cells while being on linestrings

texts:
We have a land protocol where we receive a fishnet of 1x1 km cells. Some cells are randomly choosed.
We need to put 4 points in each cell and these points have to be on a road too. The minimum distance between points have to be 500m for every points of every cells IF POSSIBLE or if it is not, we want the maximum possible distance.

In a first try we divided every cell in four 500x500 m cells with ST_CreateFishnet then we put points at the centroid of the sub-cells then on the closest road (ST_ClosestPoint). We get some good results but in the example below you can see than point 5 is too close from 6 and could be moved on the left road.

WITH
r1 AS (         -- only sub-cells which intersects random cells
    SELECT id_maille, ROW_NUMBER() OVER() AS id_grille, fishnet_500.geomgrille
    FROM fishnet_500
        JOIN t_mailles
            ON ST_Intersects(ST_Buffer(t_mailles.geom,-200), fishnet_500.geomgrille)     -- buffer < 0 to not select neightbours
)
,
r2 AS (         -- cut roads in every cells
SELECT id_maille, id_grille, ST_Intersection((ST_Dump(roads.geom)).geom, r1.geomgrille) as geomroute
FROM roads
    JOIN r1
        ON ST_Intersects(roads.geom, r1.geomgrille)
)
                -- select point on each road the closest to cell centroid
SELECT r2.id_maille, r2.id_grille, ST_ClosestPoint(ST_Union(r2.geomroute),ST_Centroid(r1.geomgrille)) as geomipa
FROM r2
    JOIN r1
        ON r2.id_grille = r1.id_grille
GROUP BY r2.id_maille, r2.id_grille, r1.geomgrille
ORDER BY r2.id_maille, r2.id_grille

If you want to give it a try I put the 3 layers (fishnet with random cells, sub-fisnet and roads) in an archive that you can find here.

I guess we can't avoid a recursive algorithm that tries many possibilities but i'm not sure.


texts: Accessing Google Tile independent from their API in OL3 is possible for paying customer?

texts:
Using map from tile providers are common use-case in web mapping: overlaying internally-hosted WMS/WFS layer on top of Google Maps/OSM/Bing Maps. These are easy in OpenLayers 2. However, this feat -especially Google Maps- is not supported anymore in OL3.
The book "OpenLayers 3 Beginner's Guide" page 107:

OpenLayers 2 had some third-party mapping APIs embedded into its core,
enabling you to use its maps inside your Openlayers-based application.
Nowadays, for better decoupling of OpenLayers API from third-party
APIs, the OpenLayers 3 team choose to have no support for other
mapping APIs that tie together tiles and an associated JavaScript API
library.
Why this decision?
One of the main goals of OpenLayers 3 was
to rewrite OpenLayers 2, making its API cleaner. The support for
Google Maps API in OpenLayers 2 has also put a significant maintenance
burden on both library and application developers, to keep up with
changes of the Google Maps API. To avoid this in OpenLayers 3, the
support for Google Maps using Google Maps API is nonexistent. However,
Google does provide its tiles independent from their API, but only to
paying customers. Fortunately, you have more alternatives to Google
Maps API (Bing, OSM).

Can anyone elaborate more on:

how OL3 can access Google Maps tile independent from their API ?
what sort of payment/contract ? Is it the same with the Google Maps API Premium ?

=================== added 6 June 2017 ==============
I found this awesome solution in OpenLayers. I tried it myself that it works with OpenLayers 4.1.1 and Leaflet 0.7.7.
Now i was wondering if this solution is the answer of my question #1 ?


texts: QGIS Visibility Analysis > Depth Below Horizon error

texts:
Background:

I need to conduct a visibility analysis that will calculate depths below the viewer's horizon.

Process:

QGIS 3.4 > Processing Toolbox > Visibility Analysis > Depth Below Horizon.

Input Data:

DEM:

32-bit floating point, single-band tif
approximately 27,000 X 17,000 cells

Observer Locations:

Point shapefile with 64 records
Geometry checks ok

Both the input DEM and shapefile are in the same CRS.

Here's a screenshot of the input:

Problem:

Running the tool results in the following error:

What does the error mean?


texts: Extract one specific value by resampling via rastercalc

texts:
I am trying to extract specific areas from the following image:

The grey areas are those of interest, they have a distinct RGB-value of 220-220-220. Thus, in my logic, the following code should leave me with the areas of interest in black and everything else white:

"image@1" != 220

Instead, I get an all black image - and I have no idea why. I tried to fiddle with brackets and conditions and stuff, but so far no difference.

/edit: Found the "solution". Yesterday, with .tiff-images, it worked. Today, with .tif-images, it does not. Good to know.


texts: How to use show Raster Attribute Table in OpenLayers?

texts:
I have some GeoTIFF files that I am showing on an openLayers map using GeoServer. Now, I wish to access a Build Raster Attribute table on my web application. In other words, if a user clicks on a point on the GeoTIFF raster layer, they should see the corresponding value shown by ArcMAP in Raster Attribute Table. I am new to the GIS.

I am thinking of creating a database with PostGIS, but I don't know how I should map the Raster Attribute Table in the database.


texts: How to include checkbox or radio button in python script for arcgis tool?

texts:
I created an arcgis custom tool for converting an excel file to both kml and shapefile .
Now I want to include radio button or checkbox to select according user interest for either kml or shapefile.


texts:
Is there any symbology set out there for sewer systems or other piped networks? The sort of elements I need to represent are:

Conduits
Junctions
Pumps
Weirs
Outfalls
Gates
etc...

There are few sets for transport systems in the public domain, but I am yet to find one for this kind of networks.


texts: Does interpolation create a vector or raster?

texts:
I have been reading about interpolation and am trying to understand whether interpolation results in a vector or raster. Say I have a layer of points and I am interested in creating a raster that fits a known grid (aligns with another raster). Do I interpolate and then convert to the raster or is raster conversion a part of the interpolation process?

Edit:

I have point vector data and a raster, so I am trying to convert my point data to a raster. From what I have read interpolation is a reasonable thing to do when converting points to a raster.

My data is stored in a PostgreSQL database, so I am trying to use Python and QGIS. I am pretty familiar with Python but not so much with this type of processing since most of the data I work with is typically tabular.


texts: Exporting from AutoCad Map 3D to postgis

texts:
Is it possible to export layers from an autocad project to a postgis database?

I'm new to autocad, so I might be mixing up my conception of a 'layer'. I have a .dwg file that contains a number of autocad layers (not display manager layers) - I understand that autocad layers do not have attributes associated them like I'm use to in GIS formats, however, each layer does have an object data table associated with it (in a 1:1 mapping). I've connected my postgis database, and I've found the 'To FDO Connection' button under the 'Output' tab, but it lists no connected data sources.


texts: Searching for imagery on Planet API: "jq command not found"

texts:
I've been trying to follow the Planet API docs to search for imagery. Specifically, when I try to make the first API call using the following code,

curl -L -H "Authorization: api-key $PL_API_KEY" \
'https://api.planet.com/data/v1/item-types' | jq '.item_types[].id'

I get the following error:

-bash: jq: command not found
  % Total    % Received % Xferd  Average Speed   Time    Time     Time  Current
                                 Dload  Upload   Total   Spent    Left  Speed
100   291  100   291    0     0    757      0 --:--:-- --:--:-- --:--:--   759
100  3034    0  3034    0     0   6552      0 --:--:-- --:--:-- --:--:--  6552
(23) Failed writing body

I understand this is probably an issue with jq, but I can't figure out why. I have double checked that all necessary packages are installed, including jq. I have also run echo $PL_API_KEY and it returns my API key as expected. Running the command with sudo does not help either. It may be worth mentioning that I installed jq with npm as opposed to pip because running pip install jq did not work for me due to build errors.
Any help?


texts:
I'm writing a little Java program to read some shapefiles & spit out the result of some calculations. I'm using GDAL, as the arcObjects.jar for ArcGIS Desktop is too restrictive for what I need (at least with my level of experience).

If I build a short Java program using GDAL imports, and I send as an .exe to my users, must my users ALSO install GDAL on their computers in order to be able to use the program? The users in question currently use ArcGIS Desktop, but not GDAL.


texts: Adding Print, Data Extract functions to published Geoprocessing widget - ArcGIS 10.1 Viewer for Flex

texts:
I have developed a pretty detailed geo-processing widget (based on a model I built and published) that takes user input, runs raster calculator and then clips the raster output to a boundary. The tool works well and the results are displayed in Flexviewer 3.5.

While I'm happy with this I now want to expand the widget functionality and include additional functions such as data extract and print functionality.  I'm wondering what is the best way of achieving this? Do I need to build this into the original model or is there a way of merging the functionality of other widgets into my geo-processing widget?


texts: Finding and merging lines that touch in PostGIS?

texts:
I would like to select and merge segments that touch like in the image below.

My table contains single segments with their own id (shown in the figure), I just want to isolate single "group" of segments that creates individual lines with their own new ID.
I used different colors to explain the different groups of segments I want to find.

I have no idea how to do that, cause I don't know the length of the lines.


texts: How to use ET GeoWizard tools or Xtools Pro in ModelBuilder or Python to automate processing?

texts: Is it possible to use ET GeoWizard or XTools Pro tools in automation of mapping in either in ModelBuilder or Python?


texts: ArcGIS Initializing License suddenly slow

texts:
I have ArcGIS 10.2.2 Desktop, License type: Advanced, installed on a windows 7 64 bit machine.

It normally takes about 12 seconds to initialize the license, whether in ArcMap or a stand alone Python script.

Suddenly today it is taking about 50 seconds. The only thing I can think of is that I ran the Advanced System Care Iobit registry cleaner but I am fairly sure I have done this before since I installed 10.2.2 without problems.

I have had similar start up times before after installing an upgraded version of ArcGIS and the only way I found to cure it was to uninstall it, remove all traces from the registry and re-install but I would dearly love to be able to fix it without going through all that again.

Can anyone suggest a way to find why it is taking so long?


texts: Measuring distance from gridsquares to contours of a shapefile

texts:
I have a baselayer (.shp file), with bathymetric contours. I want to measure the distance from each grid square in a grid, to the nearest nominated depth (so was hoping to have a way to also nomate specifically to -120 m in the .shp file).

I generally have a basic understanding of the functions of QGIS, but have absolutely no idea where to start.

I'm using QGIS v3.2.1

I have tried the below answer from Kazuhito, but appear to be getting scores that are UTM figures. I've attached photos of where I am at now.

 i.e. not a true reflection of any distance to the 100 m contour.


texts: What are the pros and cons of each mashup baselayer?

texts:
I've been making some mashups recently in OpenLayers and find myself naturally favouring OpenStreetMap as a baselayer, I think this is mostly because I like the concept. I realised the other day that I didn't really know why I prefered it, and have played about with googles various base layers, andlong with those of Yahoo and Bing.

Are there any discernible benefits of any of the baselayers that I have mentioned (or any others that I havent)? Do you find any problems with any of them or do you find that your users prefer one above the rest?


texts: Unable to get points in buffer on selected feature in QGIS using Python?

texts:
I have to create a 1.5 km buffer on each latitude/longitude of layer A. I would like to know how many points of line B will be inside the buffer of layer A. I used the following code but it takes too much time.

Is there a faster way to do this?

       crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic
       crs_bng = proj.Proj(init='epsg:4326')

       # then cast your geographic coordinate pair to the projected system

       d1=OrderedDict()
       d2=OrderedDict()
       temp=[]
       temp1=[]
       lat=0.0
       lon=0.0
       a1=0.0

       templist=[]

       Alatlong=[]
       attrs=A.getFeatures()
       hlidx = A.fieldNameIndex('Y')
       #print idx
       hloidx = A.fieldNameIndex('X')

       Alatlong=  [(i[hlidx],i[hloidx]) for i in attrs]

       Blatlong=[]
       attrs1=B.getFeatures()
       plidx = B.fieldNameIndex('Latitude')
       #print idx
       ploidx = B.fieldNameIndex('Longitude')

       Blatlong=  [(j[plidx],j[ploidx]) for j in attrs1]
       distance =1500

       for i in Alatlong:
            ans=2

            input_lat=i[0]
            input_lon=i[1]
            x, y = proj.transform(crs_wgs, crs_bng, input_lon, input_lat)

            point_1 = geometry.Point(x, y)
            circle_buffer = point_1.buffer(distance)

            for j in Blatlong:

                input_lat1=j[0]
                input_lon1= j[1]

                x1, y1 = proj.transform(crs_wgs, crs_bng, input_lon1, input_lat1)

                point_2=geometry.Point(x1,y1)

                x=[]
                if point_2.within(circle_buffer):
                   print 'True'


texts: How to EASILY add a map legend to a map in Geoserver?

texts:
I'm new to Geoserver and Neatline and somewhat to HTML coding! But, I am clueless as how to 'EASILY & Cheaply " include a simple legend on a map using a WMS statement. I have been able to use singularly the WMS "get legend request to get a simple legend, Is it possible to combined the "get legend" request with the "get map" WMS request into one WMS statement which I could imbed into a Neatline exhibit and publish the map with a legend!

Is this possible? Or if not what is the easiest way?

Thanks

Robert Jackson


texts: How to create complex shape files (polygons) in GeoJSON?

texts:
I have a project to create historical maps and would like to use GeoJSON as a default file format. However, I'm finding that there are 2 main problems with this:

There are not much available historical maps in GeoJSON/shapefiles
Creating shapes in GeoJSON is not that intuitive or accessible even with such tools as geojson.io

With regard to the latter, creating the shapes to conform to coastlines and natural boundaries is difficult. I haven't found a tool or technique (yet) that will enable me or others to produce these historical maps in GeoJSON to scale.
The strategy I am thinking of going with is to use current country Shapefiles(these resource exist, NaturalEarth) and then just modifying them to fit my needs.
Does anyone have any suggestions, thoughts, solutions on how to create historical maps in GeoJSON easily and intuitively?


texts: Territory Design based on specified criteria using ArcGIS Pro

texts:
I am looking to design territories, which take into account various criteria, such as equal population, compactness, and contiguity between basic areal units. In ArcGIS Pro, it can be achieved with Build Balanced Zones and Territory Design tools. However, they can only design territories based on a target value.
Is it possible to create territories that also take underlying boundaries into account?
Those boundaries cannot be hard criteria. They can be seen as, if possible, respect that region. I believe that it should come from a contiguity strength between basic areal units.


texts:
I am trying to draw a polyline from a single reference point. I have each specific segment of the lines Azimuth & distance. In the sequence I have areas where my line shifts over +90 / -90 of the current azimuth and then continues with the next azimuth below. I would like to preferably import these via excel and produce a record in the attribute table where I can see each segments spatial extent so I can take the Easting / Northing from the midpoint and incorporate that into an outside table. I have tried Azimuth and bearing but am having issues with it and I have to create a line for each segment. I tried advanced digitizing tools but the reference angles are from +- 90 N/S and +- 180 E/W which I could convert but seems like a hassle. Wondering if there is an easier way or if I have missed something? Here is an example database.

0.00    0.50    273
0.50    1.00    285
1.00    1.50    278
1.50    2.00    295
2.00    2.00    385 // Offset
2.00    2.30    283
2.30    2.60    283
2.60    3.05    285
3.05    3.55    285
3.55    3.85    270
3.85    4.35    270
4.35    4.75    270
4.75    5.25    270


texts: Clipping a vector line shapefile in QGIS based on a raster shapefile?

texts:
I have a vector line 'straten' shapefile which I want to have cliped based on raster shapefile "Luchtfoto 2016". Because some streets are running out of tmy raster, I want to have the vector cliped.

I know this in Arcgis but not in QGIS. How can I do this?
Yes I can easily deleted the shapes that run out of the raster, but some lines are just too long and only consists of one segment.


texts:
I'm just starting to get my head around writing simple expressions in QGIS. I have 2 polygon layers: the first shows the boundaries of houses of various types and how many people live in each. The second shows the boundaries of zones which group certain houses together.

I would like the second layer to show a summary of how many people live in each type of house in that zone. I also need the values to update automatically should i make changes to the first layer. For this i reckon that the second layer will need a different field for each house type and that these fields should be virtual.

I have gotten as far as downloading the Reffunctions plugin and getting the intersecting_geom_sum function to add up all of the people in each zone. However, i am struggling to introduce the condition whereby each field only adds the values of a specific type.


texts:
I have placed a grid onto an open source Ordnance Survey map that corresponds to Ordnance Survey Eastings and Northings. I have a value for each grid square between zero and 30.
How do I use these values to produce a density map with 4 categories in QGIS 3?


texts: QGIS Ubuntu 3.4 SQL Query Composer in WFS Connection Dialog

texts:
The QGIS 3.4 documentation mentions the procedure of filtering the features of a WFS through the connection dialog box either by double-clicking on the layer of interest, or by clicking the "Build Query" button.  However, there is no example of query building in the resulting SQL Query Composer dialog box, and I see so way to build a query in the same fashion as you would through the "Filter" function from the Layers panel after adding the WFS to the map.  Specifically, I see no way of pulling in the field values to select from.  Which might be workable if you have used the particular layer before and made note of the columns and  all their values, but unless I'm missing something, it's not possible where you are using the WFS for the first time and don't know how the database is structured.  The clickable options seem to be geared more to spatial processing functions rather than selecting field values for filtering.

What am I missing?


texts:
I have a displayed markers for the leaflet-knn on the map and the name of marker is displayed on the div at right-panel. When I click on this list, I want to show the popup on the associated marker. I have to pass selected text on right panel click but how can I show.

This is code so far I have tried.

HTML

Nearby Search Result

JS

var gjLayer = L.geoJson(testCities, {
  onEachFeature: onEachFeature
});

function onEachFeature(feature, layer) {
  layer.on('click', function(e) {
    let feature = this.feature;
    let content = "Name: " + feature.properties.name;

    e.latlng.layer.bindPopup(content);
  });
}
$(function() {
  let res;
  const distance = 10000;
  const longitude = myloc.lng;
  const latitude = myloc.lat;

  res = leafletKnn(gjLayer).nearest(
    [longitude, latitude], 5, distance);

  for (i = 0; i < res.length; i++) {
    var result = res[i].layer.feature.properties;
    map.addLayer(res[i].layer);
    searchResult(result);
  }

  function searchResult(result) {
    item = '';
    item = '
' + '
' + '' + result.name + '

';
    $(".ui.list").append(item);
  }

  $('.content').on('click', function(e) {
    let content = $(this).text();
    onEachFeature(content);
  })
});


texts: How to reflect the active state of a tool button in QGIS

texts:
I am developing a QGIS plug-in (3.0+) with a tool button (after selecting the toolbar icon you have to click in the map) and I can't set the icon active state (visual effect like the button is pushed in) like the info or the select button in the standard toolbars.
I create a new class based on QgsMapTool.
class PointTool(QgsMapTool):

    def __init__(self, canvas):
        QgsMapTool.__init__(self, canvas)
        self.canvas = canvas
        ...

    def activate(self):
        super().activate()

    def deactivate(self):
        super().deactivate()

    def canvasPressEvent(self, event):
        ...

Should I call some other method(s) in self.activate to change the active state of the icon in the toolbar?
The source code is available here: https://github.com/zsiki/arearatio


texts: Selecting features by attribute containing a string (PyQGIS)

texts:
Is it possible using Python to select all features whose attribute contain a specific string?
Since I'm a QGIS newbie I'd really appreciate your help!

The layer has an address column ("Bauvorha_3") which looks like this:

Berlin, Am Alexanderplatz, ...

Berlin, Goethestrasse, ...

Hamburg, Alsterufer, ...

Hamburg, Reeperbahn, ...

I need to select all features containing for example "Hamburg" or "ufer".

place is the user-input string.
I've tried to use LIKE but it doesn't seem to work for me.

place = self.dlg.lineEdit.text()
place = str(place)
versuch = "'{}'".format(place)
search = '"Bauvorha_3" LIKE "{}%"'.format(versuch)
self.iface.activeLayer().selectByExpression(search)


texts: What Explained Variance is good to proceed with Predictive Soil Mapping?

texts:
I have 1851 soil samples data on pH covering a study area of 7482sq.km in northern Ghana and I am using 52 environmental long-term average variables (Relief, Climate, MODIS Reflectances and derived products) to fit a model in order to explain the variability for pH prediction. So far, all models tested have shown low explained variance and sometimes even gives negative.
How may I improve the Explained Variance below?
See pictures below, the spatial distribution of points.
Below is the summary of my models explained variance.
Multiple Linear Regression: 0.03
Step-wise Multiple Linear Regression: 0.04
RandomForest: -7.02
ExtremeGradient Boosting: 0.03
Support Vector Machines with Polynomial Kernel: 0.03
Additional information about the sample raw data

Avg. distance between two sample points using nearest neighbor analysis: 2551.29m
Data source: Student research
Sampling method: Grid/Management Zone Hybrid Soil Sampling method. The grid size is 2-4sq.km, management zones are subdivisions within the grid.

My assessment so far

Removed points falling outside the study area
Removed pH value outliers
Normality check was ok (please see attached)
Variography shows a spatial structure (please see attached)
Tried recursive feature elimination to reduce the dimensionality but did not show any improvement
Tried reducing dimensionality by removing highly correlated covariates at a threshold of 0.75

I seek insights into any techniques that could help improve the model explained variance.
ATTACHED IMAGES BELOW


texts: Different number of columns and rows in raster datasets

texts:
currently I am dealing with the creation of a landuse file. I have geoprocessed a google earth .tiff-file over my DEM and created three polygon featureclasses via the editor tool containing areas of a certain landuse. After the conversion into raster files of course the number of rows and columns is different to the ones in my DEM - is there a way to fix that? (tried clipping already, no success)

In the end I'll need one single landuse file with the different categories (probably there's a easier way to achieve that!?)
I am using arcgis 10.1


texts: Is there any way to push the DPI limits when exporting a map from composer?

texts:
I noted that on OS X (El Capitan) running QGIS 2.18, I can't export map in A0 with a resolution higher than 340 dpi. If I try I receive an error message about possible memory overflows. I have 16GB of RAM on my machine so I don't think the limit is triggered by the OS.

I tested the same settings changing only the dpi. I noticed that I start getting a warning message about the memory usage at 121 dpi with a dialog box appearing just after clicking to export as jpg - I start gettting the error at 121 dpi with message To create image 7622x8829 requires about 201 MB of memory. Proceed?.

Even if I click on OK, I will then get an additional error box:

Is there any way around?


texts: Using Field Calculator to correct the nomenclature of APN values

texts:
I am currently working (in ArcMap 10.4) with two datasets that use different nomenclature for their assessor parcel numbers (APN) values. In one dataset the APN is a string dataset with the format "12345678". In the other dataset the APN is again a string but displayed as "123-456-78".

My goal is to add '-' characters to the correct position within the first string, to allow me to perform a join between the two APN fields respectively. Currently my script within the code block is seen below. I have created a new field titled "APN" and am hoping to populate this new field using the PRCL_ID string.

I am still not confident with scripting within the field calculator tool and have based this code on a similar problem in another thread. Does anyone have any suggestions on the best way to approach this problem?


texts: Where can I find the finest granularity of 2010 US census population information?

texts:
I am writing a re-districting application where I want to be able to create US voting districts based on geographic coordinate population data. In a perfect world, I would have the population for each house and the house's geographic coordinates. I am guessing that per house data is not likely to be available. What is the lowest level available? Blocks? Block groups?

I located a US census website which contains population data for block groups. Is this the finest granularity available for 2010 census data?

My end goal is to implement the shortest splitline algorithm for districting with 2010 census data. My intention is to slice up polygons using the shortest splitline according to the best population numbers relating to latitude/longitude coordinates I can find for people in the polygon.


texts: Is it possible to create a QGIS workspace with pre-loaded layers, styles and print templates with Python?

texts: I am building a python tool which aggregates data from a number of sources and generates some web-maps from this data. I was wondering if it is possible (maybe using PyQGIS) to create a QGIS3 workspace in a specified directory with the layers and styles pre-loaded?


texts: Is is possible to run the QGIS python module server side?

texts: I'm pretty new to programming but I want to use the QGIS python module on a server. Basically I am trying to create a JPEG map programatically on a server. I already have a python script to create routes server side I am just trying to find the best way to create a jpeg map with an OSM (or Google) background map, labels, and a legend. QGIS seems like it might be able to do this I am just not sure how difficult or if it's even possible to do on the server using python. I know I could probably do it via Leaflet but since I don't actually need a web map I was hoping QGIS might be a solution. Also it's a Unix server if that matters.


texts:
In my application i add layers dynamically. The first added layer is a baseLayer. When i add layer i extend bounds so the new layer could fit map (it works). The problems is that i don't know how to recalculate zoom levels. In my map config i have maxResolution set to auto and the zoomLevel is set to 0. When i add map with bigger extent the zooms are not recalculated, so i cannot see whole layer. I use GeoExt framework with OL. Code:

addLayer:function(record) {

        var layer = _factoryManager.getFactory('layer').getLayer(record, this.getServiceType());

        var map = this.getActiveMap().map;

        map.addLayer(layer);
        console.log('adding layer');
        //console.log(map.getZoomForExtent(map.maxExtent));
        //console.log(map.getZoom());
        map.maxExtent.extend(this.calculateBounds(layer));
        //console.log(map.getZoomForExtent(map.maxExtent));
        //console.log(map.getZoom());
        this.zoomToLayer(layer);
    },

And map creation:

map = Ext.create('OGC.core.panel.MapPanel',{
            map: {
                controls: controls,
                //allOverlays: true, if nothing specified all overlays is set to true
                scope:this,
                projection:this.getCRS(),
                maxExtent:bounds,
                maxResolution:"auto",
                //allOverlays: false,
                eventListeners:{
                    zoomend:function(evt){
                        console.log(this);
                        console.log('zoomend');
                        console.log(this.getZoom())
                        console.log(this.getScale());
                        Ext.getCmp('map-bottom').update('Scale: ' + this.getScale());
                    }
                }
            },
            layers: [layer]
        });


texts: Read the attributes of a WMS layer by clicking on the map

texts:
I'm trying to display a popup (containing the attribute data) by clicking on an entity in my layer with the WMSGetFeatureInfo control from OpenLayers. It's turning but without displaying the popup
Note: I'm working with OpenLayers-2.13.1 , geoext2-2.1.0,  ext-js-gpl-master
var info = Ext.create('GeoExt.Action',{
  tooltip: "Feature Info",
  icon:"Librs/geoext2-2.1.0/resources/images/default/information.png",
  text: "Get Info",
  enableToggle: true,
  allowDepress: true,
  enableToggle: true,
  control: new OpenLayers.Control.WMSGetFeatureInfo({
    //autoActivate: true,
    url: "http://localhost:8081/geoserver/STAGE/wms",
    title: 'Identify features by clicking',
    queryVisible: true,
    eventListeners: {
      getfeatureinfo: function(event) {
        if(popup) {map.removePopup(popup);}
        popup = new OpenLayers.Popup.FramedCloud(
        "getInfo",
        map.getLonLatFromPixel(event.xy),
        null,
        event.text,
        null,
        true
      );
      map.addPopup(popup, true);
      }
    }
  }),
  map: map,
  toggleGroup: 'tools'
});


texts: Extracting above ground features from unclassified LAS

texts:
What's the best way to vectorize features from LAS?

I have "Classified" LAS files, however they were never processed for classification.
Thus of 181 million points (9 tile sample) about 60% are class 1-Unassigned; 35% are class 2-Ground; and the remaining are Key-Reserved and Noise.

Is there a way to extract features (buildings, powerlines, tree canopy) using these LAS files in ArcGIS?

If this cannot be done in ArcGIS, what programs/software are optimized for such processing?


texts: Keep the name of the files in the iterations (ArcMap)

texts:
I want to use the iteration tool (ModelBuilder - insert - iterations - raster) and when I use "Zonal Geometry as Table" and write "t_%Name% some errors appear, if I use t_%n% it works but I need to keep the original name for reference the table

You can see errors here:

This one is the complete scheme:


texts: Selecting polygons based on area of overlap in QGIS?

texts:
I am working with point data in QGIS, where the points represent electricity pylons on a transmission grid.

I have created buffers with a 1 km radius around each of the points, and I would like to know if there is any way that I can extract buffers from this layer based on the area of overlap with other 1km buffers.

I only want the buffers which share less than 10% of their interior space with another 1 km buffer in the same layer.


texts:
I have a point layer with 2 columns, company name and short name. I have many company name entries but fewer short name entries with some of them being categorised under the same alias. See below for example of my data (can not show original for licensing reasons).

I have decided to go with rule-based styling because I would like to add an expression which states when "short name" = 'Ex 2' but not when "company name" = 'Example 2'. This would style Example 2 and Example 3 the same but identify Example 2 differently.

I realise I could just say:

"company name" = 'Example 3' AND "company name" = 'Example 4'

but that would take some time with my data.

EDIT

I have the answer, I can write:

"short name" = 'Ex 2' AND "company name" != 'Example 2"


texts: Sample tool in ArcGIS Spatial Analyst not working?

texts:
In ArcGIS 10.1 for Desktop I am trying to use the spatial analyst tool 'sample' to extract data from several layers.

I have double checked and the input location raster is in the same format and cell size as the input raster, however the tool will not run.

The couple of times I have got the tool to run it has not populated the table with data or the data has been incorrect.

I have not previously had a problem whilst doing something similar with other data sets on an other computer.

Any suggestions?


texts: What if format for Getcapabilities document

texts: I want to add capabilities URL for cascading a layer from remote WMS in Geoserver. What will be the exact format for Getcapabilities URL ?


texts: Assessing relative positional error in a DEM

texts:
I have a set of ground control points collected with a differential GPS and I am trying to assess the relative positional error (in contrast to absolute elevational error) in a DEM. I would like to know RMSE and standard deviation. My intuition tells me that I want to calculate all pairwise distances among ground control points and take the difference with the set of pairwise distances in the corresponding coordinates of the DEM as follows (in R):
set.seed(123)
pointNum <- 5

truth <- runif(pointNum, 70, 100) # my ground control points, elevation varying form 70 to 100m

noise <- rnorm(pointNum, 0, sd = 1) # make some noise with SD of 1

dem <- truth + noise # add the noise to simulate DEM with erros

distTruth <- as.dist(outer(truth, truth, `-`)) # calc the distance (raw difference) matrix between all control points
distDEM <- as.dist(outer(dem, dem, `-`)) # do the same for my dem

RMSE <- sqrt(mean((distTruth - distDEM)^2)) # root mean squared error between the matrices
SD <- sd(distTruth - distDEM) # standard deviation between the matrices

Is this the correct way for calculating relative positional error?


texts: QGIS processing default output vector layer extensions

texts:
You can set the default output-extension of a model in processing > options > general. The intermediate transactions in a model are processed in that default file format.

After installing QGIS, the default extension is set to 'shp'.

A well-known disadvantage of using shp is the restriction in attributenames (set to 10 characters).

When I set the default to eg. 'geojson' I can evade the naming-restriction.  But I then get some other errors on processing geometries.

So I'm looking for an explanation of the advantages and disadvantages of the different output-extensions used in the graphical modeller.


texts: QField 1.0.0: How to clip and use a Google Satellite Hybrid basemap offline?

texts:
After a few attempts I've transferred my QGIS project successfully to Qfield. I'm using Google Satellite Hybrid as my basemap and need to use it at as high a resolution as possible offline.

I've tried saving the image as a raster layer (Geotiff) in QGIS using the map view as the extent (Export - Save As - Rendered Image - Geotiff - Create VRT: off - EPSG:3857 - Extent: Map View)

To which I get an error: Unsupported Data Source: not a supported raster data source. Tried various other combinations (changing data output, VRT option) but keep getting the same error.

I've also tried to create a basemap via QField Sync. I don't get an error message but when I open the project in Qfield the basemap is in the layer list but is not displaying.

I've also tried a georeferenced image saved from Google but the resolution was too low to be usable.


texts:
S2_1 <- crop(S2_1, area)

In general it seems that this command in R works properly. But after comparing the extent of my area and S2_1- file I realised that they have not exactly the same extent.

> area@extent
class       : Extent
xmin        : 556161.6
xmax        : 604551.6
ymin        : 5305990
ymax        : 5359290

> S2_1@extent
class       : Extent
xmin        : 556160
xmax        : 604550
ymin        : 5305990
ymax        : 5359290

What can be a reason? and what can I do to avoid this? because I need to stack them :)


texts: Hyperlink relative paths to files in points layer in QGIS?

texts:
I would like to add hyperlinks to each point in a points layer in QGIS 1.7.0 running on Mac OS X 10.6.7, but I have several requirements:

(1) Each hyperlink has to be a relative path to the linked file.

(2) The linked files are mainly JPEG images, but may occasionally be HTML files. Either way, I want them to open upon click in the OS's default image/html viewer, regardless of what computer I open it on.

(3) The linked filenames are already specified in a column in the point layer's attribute table. I want to use those filenames in they hyperlinks so I don't have to manually type in a path for each point (there are several hundred points, each pointing to a different file).

How do I do this? Thanks.


texts:
i'm new to OGR processing and KML format so be gentle:

i've got a series of KML files all containing multiple regions partitioning a larger area in each.  i also have a series of geocoded points for which i want to ask WITHIN/CONTAINS queries with respect to the regions.  two questions:

i don't know how to ENUMERATE all the regions/placemarks contained in the file.  [UPDATE:] thanks to this hint i guess i can get what i want via a JSON encoding:

drv = ogr.GetDriverByName('KML')
inf = 'path/file.kml'
ds_in = drv.Open(inf)
featureTbl = {}
for kml_lyr in ds_in:
    for feat in kml_lyr:
        ftbl = eval(feat.ExportToJson())
        fname = ftbl['properties']['Name']
        geom = ftbl['geometry']
        featureTbl[fname] = geom

but there must be a better/XML-parse-like way to get "child" placemark elements?

since i know the point will only be in one of the regions, i wonder if it is worth the effort of first computing region centroids and asking the WITHIN/CONTAINS query of the regions in order of their proximity to the point, rather than just dumbly enumerating all the regions and asking until i get a hit?


texts:
I have a shapefile of contour lines, with an ALTITUDE attribute with the height ASL in metres. The source data is in 0.5m intervals, but I'd like to filter it down so I get contour lines at (for eg.) 5m intervals.

I tried doing a query on the layer using "ALTITUDE" % 5 = 0, which almost worked except I get two contour lines where I would expect one - I'm guessing this is because % is an integer operation so 245.5 % 5 and 245.0 % 5 both = 0?

Is there a better way to filter for 5m (or some other interval, say 1m or 2.5m) data?


texts: PyQGIS Problem adding certain vector layers to QGIS from postgres

texts:
I've created a tool to find layers from our PostGIS database based on a string match and add them to the current QGIS project. In the example below I'm searching for 'landscape' in the dataset name.

Problem: it adds some datasets to the current project- but not all (2 out of 4 work).

The basic code that adds the layers in: ## I feed in all the parameters with 'P' on the end

uri = QgsDataSourceUri()
uri.setConnection(hostP, "5432", databaseP, userP, passP)
uri.setDataSource (schemaP, dataChoiceNameP, geometrycolP)
vlayer=QgsVectorLayer (uri.uri(False),dataChoiceNameP,'postgres')
currentProject.addMapLayer(vlayer, False)
mygroup.addLayer(vlayer)

Note that non of these layers originally had the ID field at the beginning and even then the same 2 out of 4 worked. Then I did a drop and create (in FME using PostGIS writer) for all datasets and putting 'id' column at the front (see settings below) and still the same two datasets didn't want to load.

When adding the layer from the QGIS Browser window, PostGIS connections the two problem datasets are added successfully, also works if I use the DB manager.

All 4 datasets have the same attribute/column data types and column names, id, geom, autocad_layer, etc. So I can't see any difference in the data to explain the problem. All datasets had their geometry checked and fixed before going into the database. These are real tables, -not views- unlike a previous post.

Q1. Any ideas why those 2 datasets don't load in the python tool?

Q2. What the difference is between the QGIS GUI (Browser panel/DB Manager) add layer code and the python equivalent?


texts: How to set up wfs-t service in Geoserver using openlayers from postgis data?

texts:
I published my postgis data layers into Geoserver.
I have enabled wfs in the edit workspace > services.
Still the "layer preview" does not show any options under wfs. See the image.
How can we set up a wfst service using postgis data in Geoserver and see the map in a webpage or qgis?
Does the Geoserver include any example wms/wfs viewers within, (same as the examples folder within the openlayers) so that I can see the wms/wfs/wfs-t with multiple layer overlays? If not how can I add openlayers viewer pages in to Geoserver folder system.

(edited)
Then, where should I put the client files? Sothat I can access from http://localhost:8080/geoserver/ ?
Which directory is the web root within the Geoserver installation?


texts: Ls factor, flow accumulation, watershade, relief in Quantum GIS 2.2

texts:
I'm a beginner quantum gis. I want to make the LS factor of 30m DEM that I have merge and set the projection in wgs 84/UTM zone 49S (EPSG ; 37249). I'm having a problem when trying to create relief and watersheds in GRASS. First, I used the function " r.shaded.relief " to create relief. Secondly, I used the function " r.watershade " to create flow accumulation . But all failed . Here are the results

Any solutions?


texts: Load multi point shapefiles into postgresql

texts:
I am so curious to know how I can change the constraint to accept multipoint when I try to use pgsql to load shapefiles. I searched a lot, but I could not find the answer.

My command line is:

C:\Program Files\PostgreSQL\9.0\bin>shp2pgsql -p -c -s 2400 -W LATIN1 C:\Projects\Mandana\OBO\Oversiktskartans_texter.shp | psql -U postgres -d postgis

I would appreciate if you tell me how we can add constraint here to accept multipoint.


texts: OpenLayers GeoJSON: limit significant digits to 6 when using write()?

texts:
Is there a way that OpenLayers.Format.JSON write() method to serialize coordinates using 6 decimals only ?

Example: instead of http://pastie.org/4200997 to have http://pastie.org/4200994


texts: How to overcome MemoryError from numpy.unique

texts:
I am using Numpy in version 1.11.1 and have to deal with an two-dimensional array of

my_arr.shape = (25000, 25000)

All values are integer, and I need a unique list of the arrays values. When using lst = np.unique(my_arr) I am getting:

Traceback (most recent call last):
  File "<pyshell#38>", line 1, in <module>
    palette = np.unique(arr)
  File "c:\Python27\lib\site-packages\numpy\lib\arraysetops.py", line 176, in unique
    ar = np.asanyarray(ar).flatten()
MemoryError

My machine has only 8 GB RAM, but I tried it with another machine with 16 GB RAM, and the result is the same. Monitoring the memory and CPU usage doesn't show that the problems are related to RAM or CPU.


texts: What lat/long formats can I use for a map data file

texts:
I can't find documentation as to the formats allowed for lat/long in data files for cartoDB maps? Should Lat/long be in decimal form, as deg/min/sec or geocoded as lat/long?

IOW is there a list of lat/long formats that cartoDB will understand if passed in a CSV or JSON file?

I think there's a presumption that is one is starting from a GIS context and thus already know the answer. I'm not!


texts:
i have around thousands of satellite images and i want to make an application for the quick visualization of these image some what like Google earth i.e the images changes when the extent is changes (Zoom in = High resolution, Zoom out = Low resolution imagery) something like that.
I have server with xeon hexa core processors with 128 GB Ram but i am bit confuse about graphics card, I have 2 GB dedicated Graphic card for this, is this sufficient. or should i go for higher for images rendering?
Kindly suggest me something.


texts: Assigning bars to closest building using PostGIS?

texts:
I'm using PostGIS and OpenStreetMap, I tried to assign some bars to the nearest building from themselves I based myself on an example, see: link You can see below picture that resulted from below script. Where its circled in red the location of each bar is incorrect. Its not showing what I need to see. The selected location of the bar is not in front of it.

SELECT DISTINCT ON (b.osm_id) b.osm_id as osm_bar, a.osm_id as osm_building, ST_Distance(a.way, b.way)  as distance
FROM buildings a, bar_not_assigned b
WHERE ST_DWithin(a.way, b.way, 10)
ORDER BY b.osm_id, a.osm_id, ST_Distance(a.way, b.way);

When i used the script in this post link
 I got this result:

select poly.osm_id
from (select osm_id, way from buildings) poly
cross join lateral
  (select osm_id, way
   from bar_not_assigned pt
   where st_dwithin(pt.way, poly.way, 10)
   order by pt.way <-> poly.way limit 1) pts;

I changed limit to 1 to have one bulding for each bar, but does not work.


texts:
I am trying to run a simple python script that changes a feature to a polygon.

However, I want the outputted polygon to have the same name an the inputted feature. The code that I ran, outputs the name as "lyr", when I really want the name to be the input name. I used "lyr" as a variable name in for the for loop

This is what I have:

import arcpy
arcpy.env.workspace = "C:\Users\Daimon Nurse\Desktop\Grounds Project"
mxd = arcpy.mapping.MapDocument(r"C:\Users\Daimon Nurse\Desktop\Grounds Project\ZonesMAP.mxd")
Layerlist = arcpy.mapping.ListLayers(mxd)

for lyr in arcpy.mapping.ListLayers(mxd):
    output = r'C:\Users\Daimon Nurse\Desktop\Grounds Project\DFMGROUNDS.gdb'
    arcpy.FeatureToPolygon_management(lyr, output + '\lyr', '#', 'ATTRIBUTES', '#')
    del lyr

del mxd


texts: Creating polygons from points using lines as barriers

texts:
This question is similar to a previous post where the member asked if it were possible to create Theissen (Voronoi) polygons from lines:

Creating Thiessen (Voronoi) polygons using lines (rather than points) as the input features?

The issue I'm trying to resolve is that I have hundreds of thousands of point features that belong to an assigned territory by clerks who may or may not be familiar with the geography of these areas. As such, many points get assigned the wrong territory value. In order to fix this, I want to find points that don't belong within a cluster and generate polygons from the points that are aligned with streets. If a street dead ends then I want the voronoi polygon to end when it meets the edge of the closest point's voronoi polygon. Once I am confident that the territory boundaries have been accurately defined, I would like to perform a spatial join and carry over the polygon territory values to the point features.

I grasp that the best approach is to convert the polygons and line features to raster in order to increase precision of Voronoi cell boundaries (although I could be totally wrong on this).

The only difference is that I want to create Voronoi polygons from point features but want the boundaries to stop when they encounter the boundary generated from another point within the feature class as well as a line feature.

This is what I want to fix so that the polygon boundaries do not cross the line:

I used QGIS to create the Voronoi polygons because I don't have the license level to create Thiessen polygons in ArcGIS. I haven't seen an option to use a barrier in create Voronoi polygons in QGIS.

I brought the features in from QGIS to ArcMap to use features to raster tool and got this:

I'm not sure where to go from here.


texts:
I'm using openlayers 3 to create web application that allow to user to draw LineString features in the map.
this is the code:

var raster = new ol.layer.Tile({
source: new ol.source.MapQuest({ layer: 'sat' })
});

var source = new ol.source.Vector();

var vector = new ol.layer.Vector({
    name: 'my_vectorlayer',
    source: source,
    style: new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: '#ffcc33',
            width: 5
        })
    })
});

var map = new ol.Map({
    layers: [raster, vector],
    target: document.getElementById('map'),
    view: new ol.View2D({
        center: [-11000000, 4600000],
        zoom: 4
    })

});
var draw;
function addInteraction() {
 map.removeInteraction(singleClick);
    draw = new ol.interaction.Draw({
        source: source,
        type: ("LineString")
    });
    map.addInteraction(draw);
}

by the previous code i can drow lines in to the map. the drawn lines will be added to  vector layer. I wont when the user select one of the line that he draw can remove them.
this is the code of selecting the feature:

var singleClick = new ol.interaction.Select();
function addSelect() {
    map.removeInteraction(draw);
    map.addInteraction(singleClick);
}

and it's work very will,

I just want the user to be able to remove the selected LineString...